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

The accurate measurement of ecosystem biomass is of great importance in scientific, resource management and energy sectors. In particular, biomass is a direct measurement of carbon storage within an ecosystem and of great importance for carbon cycle science and carbon emission mitigation. Remote Sensing is the most accurate tool for global biomass measurements because of the ability to measure large areas. Current biomass estimates are derived primarily from ground-based samples, as compiled and reported in inventories and ecosystem samples. By using remote sensing technologies, we are able to scale up the sample values and supply wall to wall mapping of biomass. Three separate remote sensing technologies are available today to measure ecosystem biomass: passive optical, radar, and lidar. There are many measurement methodologies that range from the application driven to the most technologically cutting-edge. The goal of this book is to address the newest developments in biomass measurements, sensor development, field measurements and modeling. The chapters in this book are separated into five main sections.

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 CO 2 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 (Borsboom et al., 2002;Domínguez, 2002;IDAE, 2005aIDAE, , 2007. 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, 2005bIDAE, , 2007Velá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 the 2010 Plan (MITYC, 2010). 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 complete review of the recent advances and future developments in this issue. 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).

www.intechopen.com
Thus, remotely sensed data have not only facilitated an increase in the speed, cost efficiency, precision, and timeliness of inventories, but they have also allowed the construction of maps of forest attributes with spatial resolutions and accuracies that were not feasible even a few years ago (McRoberts & Tomppo, 2007). These advances have made remotely sensed data the primary source for biomass estimation (Lu, 2006). However, no studies have yet presented a technique for biomass estimation that is consistent, reproducible and entirely applicable at regional scales (Muukkonen and Heiskanen, 2005;Powell et al., 2010), especially in forests located in areas of irregular topography and characterized by heterogeneity in species composition, complex stand structures and environmental conditions (Brown et al., 1999;Lu & Batistiella, 2005;Lu, 2006;Mallinis et al., 2004;). Those characteristics are inherent features of Mediterranean forests (Shoshany, 2000), where not many studies have been carried out. The objective of this chapter is to explain a methodology developed to estimate the amount of FRB potentially suitable for renewable energy production in the pine forests of Mediterranean areas at regional scale, using satellite images and forest inventory data. It is intended, therefore, to eliminate a major barrier to the use of this renewable source of energy. In turn, by using a plain methodology, it is intended that the method developed can be adopted by decision makers and land managers for both forest management and regional planning, considering that energy planning is a major component of land management. For this study, we used FRB data obtained from allometric equations applied to the Second Spanish National Forest Inventory (NFI-2) (dependent variable) and spectral data from Landsat 5 TM imagery (independent variable). In order to avoid the effects of forest heterogeneity in the establishment of accurate predictive models, different methods were tested to extract the spectral data from the Landsat images.

Materials and methods
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.

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 www.intechopen.com (2002)(2003)(2004)(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 km 2 ), (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).

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.

www.intechopen.com
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) ( 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.

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., 2001Foody et al., , 2003Gasparri et al., 2010;Hall et al., 2006;Labrecque et al., 2003Labrecque et al., , 2006Lu, 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, 1998bSteininger, 2000;Tangki & Chappell, 2008;Wulder et al., 2008;Zheng et al., 2004Zheng et al., , 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 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.

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 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, 1998bShoshany, 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.

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;andBatistiella, 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): 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.

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 www.intechopen.com 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).

Spectral segmentation
The main problems in constructing predictive models related to the extraction of radiometric data using high-resolution aerial photographs of delimited forest area are: (i) human errors in the establishment of limits during the visual image interpretation process, (ii) assumption of a FRB constant value for the entire homogeneous area based on a point location that is not representative of larger forest stands (Pekkarinen, 2002). Image segmentation can be applied to delimit homogeneous spectral features as a reference for spectral data extraction (Hall et al., 2006;Mäkelä & Pekkarinen, 2001;Pekkarinen, 2002). In this study, the segmentation algorithm RGB clustering incorporated in the software Erdas Imagine was applied. The essential parameters that control the results in this method are the RGB composition (which is used to apply the segmentation), the stretch applied to each band in the composition, and the number of bins into which each band is divided. Four different segmentations were performed, modifying the stretch (2 or 4 standard deviations) and the bins (7-6-6 and 4-3-3) using a RGB composition TM7-TM4-TM3 to model the heterogeneity of the analyzed forest (Table 2). In addition, in order to remove pixels within the immediate vicinity that were not attributable to the IFN-2 plot, the obtained image segments containing field plots were then intersected with 3 x 3 pixel windows, resulting in a mask of homogeneous pixels. Thus, the mean pixel value for each image band was then computed from the pixels that belong to the same spectral category as that of the central pixel of the NFI-2 plot ( Figure 5).

FRB mapping and validation method
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 respective estimation model and its validation each time.
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 (R 2 ) was used, and the Root Mean Square Error (RMSE) and the relative RMSE (RMSE r ) 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 RMSE r were calculated to evaluate the results.

Model application and inventory
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 R 2 and RMSE and RMSE r , 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 resource and its spatial distribution in the study area. 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 and a decreasing number of plots. All of the groups showed significant correlations, with the majority of the considered variables yielding p-values of generally p < 0.05, with the one exception being the first CV group, as this group contained only three plots. In all groups, the highest coefficients of correlation were obtained for variables related to wetness information (TM band 5 -TM5-, TM band 7 -TM7-, the third principal component -PC3-, the third tasseled cap component -TC3-, the moisture stress index -MSI-, and the sum of TM5 and TM7 -MID57-), reaching similar coefficients than the Normalized Difference Vegetation Index (NDVI). The coefficients of these variables increased from values of 0.450-0.460 in the group of percentile 10 to more than 0.850 in the group of percentile 3 with all results being statistically significant. These results show that the degree of spectral heterogeneity determines the feasibility of building accurate predictive models. Thus, if the groups with more plots are used, the models will have low prediction capacity. By contrast, if the groups with lower numbers of plots are used, the models will have higher predictive capacity, but, in turn, the models will be biased to the sample, being not representative of all of the forests of Teruel Province. Therefore, we selected the CV4 group because, among the groups characterized by an elevated homogeneity of plots (allowing execution of models with high R 2 ), it is the one with the most number of plots (68). Scatter plots comparing FRB and the spectral variables clearly reveal nonlinear relationships. As a result, the most suitable standard transformation was applied to the independent and dependent variables in order to guarantee the linearity principle in the multiple linear regression model (logarithmic -ln-, square root -sq-, exponential -exp-and inverse -inv-) (Hair et al., 1998). Table 4 shows linear regression models performed with each one of the five random divisions of its sample. All of them include only one variable, which is always related to water content in vegetation. The lack of a second variable in these models is due to the high degree of auto-correlations between the spectral variables. This avoids violation of one of the principles of multiple linear regression models. Even so, the R 2 values achieved were higher than 0.7. Consequently, they are suitable for FRB estimation in our study area. The results obtained from validation using 20% of the sample showed that the model run with MID57 (TM5+TM7) is the most suitable, because the RMSE r was only 26.67%, significantly lower than the others. As a result, it was selected to obtain FRB cartography with the 1994 Landsat image, using the pine forest areas in the Aragon 1:50,000 cartography as a mask (Figure 6).  Table 4. Linear regression models obtained from 3 x 3 pixel windows Fig. 6. The regression model selected using fixed pixel windows and the FRB map produced applying it

Models completed using visual analysis to extract the radiometric information
Correlation analysis completed using homogeneous areas derived from visual analysis showed that spectral variables related to wetness yielded the highest Pearson coefficients. Among them, TM5, TM7, TC3, MSI and MID57 were identified as the best predictor variables, with an R > 0.68. In addition, vegetation indices were again the second group of variables that showed high correlation levels, with NDVI showing the strongest correlation (Table 5). Finally, note that scatter plots comparing FRB and spectral variables also showed nonlinear relationships; therefore both independent and dependent variables were transformed using standard transformation before introducing them in the SPSS software for the linear regression analysis.  Table 5. Pearson's coefficient of correlation (R) between spectral variables and homogeneous areas delimited by visual analysis (* correlation is significant at the 0.05 level)

Variable
Consequently, it is possible to run estimation models using visual homogeneous areas with similar accuracy level as with fixed windows because the correlation coefficients are comparable. Therefore, the main difference in the analysis made using visual homogeneous data comes from the number of plots included, which is nearly double that of the group of percentile 4 (131 versus 68). As a result, the probability of constructing an over-fitted model is reduced (Hair et al., 1998). Therefore, this method to extract radiometric data appears to be more suitable to perform regression models for FRB estimation, as these models could be more representative of all environmental characteristics of Teruel pine forests. Only one of the five models successfully included more than one independent variable (Table 6), owing to the high auto-correlation between them. The variable chosen in each regression attempt depends on sample divisions; but in all cases, this variable was related to wetness. The model run with the sample N 3 was selected to produce a map (Figure 7), as it showed good conciliation between its R 2 and its RMSE r . In addition, this model allowed direct comparison to the previous extraction method because both models use the same radiometric variable (MID57).  Table 6. Linear regression models obtained using homogeneous areas delimited by visual analysis Fig. 7. The regression model selected using homogeneous areas delimited by visual analysis and the FRB map produced applying it

Models carried out using spectral segmentation to extract the radiometric information
The correlation analysis using the data derived from the third extraction method was only performed using the plot groups CV7 and CV6 for two reasons: (i) to directly reject in the analysis the plots with high probability of containing radiometric data related to different spectral features than the FRB data; and (ii) to identify a regression model with similar R 2 values than were derived from the use of fixed 3 x 3 windows, but using a higher number of plots, making it more representative of the entire study area. In addition, only NFI-2 plots inside areas classified as pine forest in the Aragon 1:50,000 forest map were considered for this analysis. As is shown in Table 7, the radiometric variables related to wetness were again the most correlated. In addition, it was shown that the correlation coefficients increase as the homogeneity in the sample increases, independently of the performed segmentation. Focusing on these results, only the S4 segmentation in group CV6 shows higher regression coefficients than using a fixed 3x3 window without restrictions. Concretely, the maximum R value S4 was reached using the variables TM7 and MID57 (R>0.7). Nonlinear relationships between the dependent variable (FRB) and the independent variables (radiometric data) were again revealed; thus transformations were applied. Table 8 shows the regression models that were obtained using plots within group CV6. Four models were performed using a wetness radiometric variable. This type of variable was also selected initially in the only one model that used more than one variable. The R 2 coefficients were situated between 0.53 and 0.59 and only one model indicated an RMSE r lower than 40% in its validation. Owing to the limited difference in terms of R 2 , this model with the lowest RMSE r was selected to derive FRB cartography for the study area ( Figure 8).

Segmentation S1
Segmentation S2   Table 7. Pearson's coefficient of correlation (R) between spectral variables and data obtained using segmentation and 3x3 pixel windows with restrictions and without restrictions (* correlation is significant at the 0.05 level)

FRB cartography validation
The accuracy assessments of every regression equation (RMSE and RMSE r ) were completed using plots that showed the same homogeneity criteria as was used to run the model. However, since the selected estimation models have been applied to each one of the Landsat pixels located in forested areas in Teruel Province, the degree of success in these must also be evaluated at that scale.
To accomplish this, the NFI-2 plots excluded from the estimation models and their validation were considered. In order to guarantee the results, those plots that were affected by inaccuracies in their field location and/or by the radiometric response of different landscape elements located in their immediate vicinity, were removed from the validation sample. Consequently, group CV8 was used since it includes a high number of plots, which ensures that the validation results were not biased by using only the ideal plots.
As it can be seen in figure 9, the results show few differences between the three maps. Those obtained from 3 x 3 fixed windows yieded a RMSE r of 64.26%, while spectral homogeneous forest areas had RMSE r values of 66, 71% and 65.06%, respectively. These results at pixel level can be considered tolerable for the study area considering previous experiments using similar methodologies for boreal environments less affected by heterogeneity than Mediterranean forests. Thus, Tokola et al. (1996), Tokola & Heikkila (1997), Mäkkelä & Pekkarien (2001) and Katila & Tomppo (2001) reported RMSE r to estimate forest parameters such us timber volume or total volume from about 65% to more than 100%. In this respect, it is important to emphasize that estimation error in cartography derived from satellite images decreases with an increase in the size of the area used to validate it. For example, Fazakas et al. (1999) showed a RMSE r of 66.5% at pixel level, but when using an aggregation area of 598 ha, the RMSE r was 8.7%. However, it was not possible to carry out a similar analysis in our study area because no other FRB data were available at any scale. Fig. 9. RMSE r differences between FRB maps obtained using the three extraction methods

Inventory of potentially available FRB
Previous research has shown the utility of wetness variables obtained from Landsat TM images to estimate FRB, regardless of the image date in the summer period. This is because these kinds of variables yield high and similar estimation FRB models using June, July and August scenes. In addition, statistical differences were not found in the moisture content of the four pines studied over the 3 month period (García-Martín et al., 2008a). In relation to results in previous sections, the model selected using the extraction of 3 x 3 fixed windows was applied to the MID57 neocanal image derived from the July 2008 Landsat image in order to inventory the FRB. This model was finally selected because it had higher predictive power (R 2 of 0.711 versus 0.595 on the model selected in the visual analysis method and 0.535 in the segmentation method) and allows the development of maps with the lowest estimation error. This last point shows that the limited set of plots used in the 3 x 3 fixed windows method is as representative for FRB estimation as those used in the other two methods, although they were composed of nearly double the number of plots. Figure 10 shows the estimation cartography obtained for the entire study area. This cartography allows calculation of the total amount of FRB resource at the provincial level, which amounts to 5,449,252 tons. In addition, with the high spatial resolution (25 x 25 m), the cartography precisely reveals the richest regions and FRB distribution within them. This makes this cartography especially suitable for determining optimal areas, taking into account other spatial variables that also determine the technical and economic feasibility in the harvest of this renewable energy resource, for example: (i) slope, which influences the possibility of using machinery and its efficiency; (ii) distance to forest tracks, which determines a portion of the transport costs; and (iii) area of forest stands, which is related to the necessary displacement during the working day. These three spatial factors and the quantity of FRB derived from remote sensing data at 25 m resolution can be integrated into a Geographical Information System (GIS) to identify areas more suitable for harvest of FRB, with attention to principles of sustainable ecological forest management (Pascual et al., 2007;García-Martín et al., 2011).

Conclusions
This study demonstrates the utility of Landsat TM images and forest inventory data in estimating FRB in Mediterranean areas. The methodology employed provides a continuous and complete estimation of FRB that overcomes problems associated with the limitations of point inventory. This information can be very useful in determining the most suitable areas in which to install power stations that make use of this resource; the lack of this type of information is one of the main problems currently facing the industry. As a result, an increase in the use of this kind of biomass resource will help to achieve the stated objectives of renewable energy production in Spain. Moreover, this is especially important considering two facts: (i) the current socio-economic context, with an increase in the price of petrol due to international instability, and public concerns regarding nuclear power, owing to Fukushima incident; and (ii) biomass is currently the only renewable energy that can be used as a strategic source of energy, as it is always available independent of weather conditions and can be stored easily prior to use (Jarabo, 1999). The three sources of information typically used for biomass estimation (data from field sampling, satellite imagery and ancillary data) (Lu, 2006) have been carefully integrated, bearing in mind the specific characteristics of the Mediterranean environment. Thus, to reduce the problem that heterogeneity introduced into estimation models, three different methods were tested to extract radiometric data from a TM image. In agreement with previous studies focused on AGB estimation (i.e. Foody, 2001;Lu et al., 2004;Lu, 2005;Lu & Batistiella, 2005;Steininger, 2000) and with previous analyses carried out with our data (García-Martín et al., 2006, 2008a, 2008b, the derived spectral variables related to wetness showed the strongest correlations to FRB, independent of the method used to extract the radiometric data. Analyses undertaken revealed statistically significant correlations with the three methods, but the use of 3 x 3 pixel windows and CV were indicated as the most appropriate to isolate a homogeneous group of plots that allow the most accurate regression equations (R 2 >0.7). In addition, despite the fact that these models were performed using half of the plots than those derived from the other two extraction methods, it was representative of the entire study area, as was shown in the validation results obtained at pixel level (RMSE r was similar in the cartographies obtained from the three methods). The visual analysis method did not yield such good results probably because of inaccuracies in the delineation of the homogeneous areas and because of the assumption of a constant value of BRF for the entire area from one precise location. The segmentation method did not improve the results either. This fact can be related to the lack of success to model effectively the spatial distribution pattern of BRF in the study area using RGB clustering. Finally, it is important to note that any of the methodologies tested to extract radiometric data successfully completed multivariate regressions models. This was because of the high autocorrelation among the well-correlated radiometric variables with FRB (wetness variables and vegetation indices). The use of these variables together will run a model affected by collinearity.
Validation cartography was based at the pixel level using plots that had not been previously used in the regression equation, removing only those with higher probability of error. This approach ensured the independence of the validation. RMSE r obtained in the three cartographies (above 65%) were better or at the same level than those produced in previous experiments conducted in boreal environments to estimate forest parameters using Landsat images. This is a positive result considering the higher complexity of Mediterranean forests. The remaining errors in the estimates of FRB can be related to the factors highlighted in García-Martín et al. (2008a,b): (i) inaccuracies in the fieldwork undertaken to establish the allometric equations and statistical analyses of the equations; (ii) problems involved in relating NFI plots to satellite data, mainly related to inaccuracies in the plot field placement (the three methodologies applied in the present study in linking ground data and remotely sensed data helped to reduce this problem, but it is still possible that errors were present in the final sample data); (iii) inaccuracies related to heterogeneity of the sample, because different pine species were considered and they are distributed in different regions of the study area; and (iv) limitations related to the spectral, radiometric, and spatial resolution of the TM in highly heterogeneous environments. However, despite these limitations, the size of the scenes and ease of distribution at no cost, make Landsat TM the most suitable in terms of achieving the objective of developing a useful methodology for estimations of FRB at regional-scales. Finally, in order to improve outcomes, different lines of research should be considered. Firstly, different segmentation methods other than RGB clustering should be explored to determine if they can better model the forest spatial pattern and, as a result, obtain more accurate estimation models. Concerning this, the eCognition and the Definiens Developer segmentation procedures offer the possibility of considering additional features together with spectral data, based on similarities in shape and size. Secondly, focusing only on Landsat data, it is necessary to explore the use of other statistical methods that allow highly auto-correlated dependent variables to be considered jointly. In addition, it would be useful to explore the capability of hyperespectral sensors to identify narrow spectral bands highly correlated to FRB and poorly correlated between them and to examine the capability of SAR data from different sensors. This will allow the execution of better predictive multivariate regression models not affected by collinearity. Lastly, it would be useful to integrate data from physical variables that can be related to FRB quantity and distribution such as elevation, slope and aspect, as well as diverse biophysical parameters such as soil, lithology and climate.