Assessment of biomass and carbon content in a Mediterranean Aleppo pine forest using ALS data

Tree biomass estimate is essential for carbon accounting, bioenergy feasibility studies and forest sustainable management. However, little research has been conducted focusing on the use of Airborne Laser Scanning (ALS) technology in Mediterranean Aleppo pine forest. Thus, the availability of ALS information provided by the Spanish National Plan for Aerial Orthophotography (PNOA) determined the main objective of this research. It is our aim to test the suitability of the low point density, discrete, multiplereturn, ALS data, to estimate and map the total biomass (TB) and its carbon content in Pinus halepensis Mill. forest. Accordingly, TB was calculated in 45 field plots, located north-eastern Spain, using allometric equations. This information was related through a multivariate linear regression analysis with a collection of independent variables extracted from the ALS data. The predictive model was validated using a leave-one-out crossvalidation (LOOCV) technique. Then, a regular grid with cell size 25 m x 25 m was generated in order to compute TB at stand level. Afterwards, biomass was transformed to carbon content by using a conversion factor. Maximum height, kurtosis and percentage of returns above 1 m were the ALS metrics included in the fitted model, which presented a R 2 value of 0.89. OPEN ACCESS


Introduction
In the last decades, operational collection of information on relevant characteristics of forest ecosystems by means of pure ground-based field inventories has been revolutionized by the development of remote sensing sensors [1].Optical and radar remote sensing have been widely used to map forest structural attributes and biophysical parameters [2][3][4][5].In this sense, Light Detection and Ranging (LiDAR) has emerged as a promising techniques in forest attribute estimation [1,6,7].Most commercial LiDAR systems are small-footprint, discrete-return airborne lasers, also referred to as Airborne Laser Scanning (ALS).These systems are able to accurately characterize the threedimensional structure of the forest canopy due to its capacity to record the different reflections of each emitted pulse, corresponding to both, vegetation and terrain beneath it [6,8].Accordingly, in order to generate a precise Digital Elevation Model (DEM), it is necessary to filter the ALS data and to interpolate the points classified as ground.These DEMs are needed to normalize the heights of the ALS point cloud with respect to the terrain heights [9][10][11].
Particularly, forest biomass mapping has gained interest in recent years due to the following three reasons: i) the relevant role of forests in the carbon cycle and in the balance of greenhouse gases emissions, ii) the assessment of the available resources for bioenergy production and iii) the sustainable management of forests [25].The total biomass (TB) refers to the dry weight of the plant material from trees, including roots, stems, bark, branches and leaves from the ground to the apex [1].Conventional methods used to estimate this variable are based on field measurements.However, these methods are limited by the high cost of such destructive sampling, necessary to create allometric models required to extrapolate the results.Alternatively, remote sensing has allowed forest biomass mapping in a wide range of spatial and temporal scales, reducing costs and fieldwork.Numerous studies have correlated the biomass and the spectral response of vegetation using passive optical sensors.However, the results have been typically affected by saturation problems when the biomass is greater than 100 tons/ha [25].On the other hand, approaches using active sensors such as Synthetic Aperture Radar (SAR) have shown that they are more sensitive to higher levels of biomass and, more recently, the use of ALS technology is becoming a tool with great potential to estimate biomass [26,1].This fact, together with the availability of ALS data captured for the entire Spanish territory under the National Plan for Aerial Orthophotography (PNOA), has determined the aim of this work, which is to evaluate the suitability of this information to estimate and map TB and carbon content in Pinus halepensis Mill.forest at stand level.Our objective implies developing a statistical model that adequately relates the information provided by the ALS data with the one extracted directly from the field.

Study area
The study area corresponds to monospecific stands of Pinus halepensis Mill.and it is located in the central Ebro valley (41° 50' N, 0º 57' W) (Fig. 2).This river crosses the Autonomous Region of Aragón, sited northeastern Spain.The forest under study is fragmented in stands of variable sizes and occupies 8266 ha.In some areas, Aleppo pine forest is interspersed with evergreen shrubs, dominated by Quercus coccifera L., Juniperus oxycedrus L. subsp.macrocarpa (Sibth.& Sm.) Ball, and Thymnus vulgaris L. Part of the study area is located inside the Military Training Center (CENAD) "San Gregorio", involving a direct risk of fire.According to the Third Spanish National Forest Inventory (NFI 3), pine forest is the most significant habitat in Autonomous Region of Aragón.In fact, it represents 49.88% of the forested area.In addition, Aleppo pines play an important role in the protection and restoration of forest as this tree species is practically the only one adapted to the adverse climatic and edaphic conditions of the study area.In this regard, it is important to notice that the study area is characterized by nutrient-poor, gypsiferous soils.In addition, this area presents a hilly topography, with elevations ranging from about 400 m to 750 m above sea level.Climate of the region is Mediterranean with continental features, i.e., irregular annual precipitation, cold winters and hot, dry summers [27].

ALS data
The ALS data were provided by the Spanish National Plan for Aerial Orthophotography (PNOA, http://www.ign.es/PNOA/vuelo_lidar.html).The information was captured in several surveys conducted between January and February 2011.A Leica ALS60 sensor was used.This is a smallfootprint, discrete return sensor.Data were delivered in 2 km x 2 km tiles of unclassified points in LAS binary file, format v. 1.2, containing x, y, z coordinates (UTM Zone 30 ETRS 1989), with up to four returns recorded per pulse.The flying height of the ALS mission was around 3,000 m.a.g.l.The sensor was operating in 1.064 μm wavelength, 0.22 mrad beam divergence and ±29 scan angle degrees.The resulting ALS point density for the study area was 1 point/m 2 with a vertical accuracy higher than 0.20 m.
After eliminating the noise of the point cloud, it was necessary to perform a filtering process to identify ground returns.The classification algorithm implemented in the MCC software 2.1, developed by Evans and Hudak [28], was applied.This algorithm was selected given its good performance in hilly forest landscapes [9].Points classified as ground were interpolated using the "Point-TIN-Raster" interpolation method [29] in ArcGIS 10.2 software (ESRI, Redlands, CA, USA) to create a DEM.This model, with 1 m spatial resolution, was used to normalize the heights of the point cloud.Then, using the "ClipData" and "CloudMetrics" commands implemented in FUSION LDV 3.30 software [30], several ALS point cloud metrics commonly used in vegetation structure modeling were calculated at each plot.A threshold of 1 m above the ground was applied in the selection of the points used in the calculation in order to exclude from the metrics the returns belonging to bare-earth and understory [31].

Field data
The fieldwork campaign was conducted between July and September 2014.The location of 45 circular plots with 15 m radius was selected via a stratified random sampling, trying to obtain a representative sample of the variability of forest height and topographic complexity [31].The selected sites were staked out in field using a Leica VIVA GS15 CS10 real-time kinematic Global Navigation Satellite System, achieving an average planimetric accuracy of 0.15 m.The total tree height (h) was measured using a Vertex instrument for precision height (Haglöf Sweden®).In addition, the tree diameter at breast height (dbh), considering the standard height of 1.3 m used in Europe, was measured using a calliper Haglöf Mantax Precision Blue (Haglöf Sweden®).In total 2,063 trees were inventoried (only trees greater than 7.5 cm dbh).The different biomass fractions at each plot was calculated using the Pinus halepensis Mill allometric models reported by Ruiz-Peinado et al. [32] (equations 1, 2, 3, 4 and 5): W b2+n = 6.197 + 0.00932 Where Ws is the biomass weight of the stem fraction (kg), Wb7 is the biomass weight of the thick branch fraction (diameter larger than 7 cm) (kg), W b2-7 is the biomass weight of medium branch fraction (diameter between 2 and 7 cm) (kg), W b2+n is the biomass weight of thin branch fraction (diameter smaller than 2 cm) with needles (kg), and W r is the biomass weight of the belowground fraction (kg).
The results of these equations were aggregated to calculate the TB in kilograms of each tree in the field plot, in order to extrapolate the TB values from sample plots to per hectare biomass value (kg of dry biomass per ha).

Predictive model for estimating the TB
Following Means et al. [12], Naesset [13], Gonzalez-Ferreiro et al. [11] and García et al. [10], a multivariate linear regression approach was adopted in order to develop a model able to estimate the TB.The statistical analysis was performed in the R environment (http://www.r-project.org/).Given the large number of potential ALS-derived metrics, the Spearman's rank correlation coefficient (Rho) was previously applied in order to select those independent variables with the strongest correlation coefficient with TB [22].Then, the selected variables were included in a stepwise regression, trying to develop a parsimonious model in order to avoid over-fitting [33,34].Predictor variables with a significance value of partial F statistic greater than 0.05 were removed from the model [31].Besides, the fitted model was selected in compliance with the statistical assumptions of linearity, normality of the residuals, homoscedasticity, and independence or no autocorrelation [10].
A leave-one-out cross-validation (LOOCV) technique was applied with the purpose of performing an unbiased assessment of the predictive capacity of the model [14].The model goodness-of-fit was assessed comparing the adjusted coefficient of determination (R 2 ), the root-mean-square error (RMSE) and the mean of the residuals obtained from LOOCV technique [35].The final model was obtained by averaging the model coefficients computed at each step of the cross-validation performed [14].

Mapping TB and carbon content
The pixel size selected to compute the ALS-derived metrics and to map TB was 25 m x 25 m, representing an area of 625 m 2 , similar to the field plot dimensions (706.86 m 2 ).Metrics involved in the equation of the TB model were obtained using the "GridMetrics" and "CSV2Grid" commands implemented in FUSION LDV 3.30 [30].Then, using map algebra in ArcGIS 10.2 (ESRI, Redlands, CA, USA), model coefficients were applied to the raster layers to achieve a TB map at stand level.Finally, TB values were converted to carbon content using a conversion factor of 0.499 following Montero et al. [36].

Results and Discussion
A summary of the variables estimated in the field is presented in Table 1.Terrain slope in the plots averaged 10.7° and ranged from 0.7 to 25.7°.As can be observed, structural variables related to Aleppo pine forest show a remarkable variability.Inventoried trees range from smaller trees (heights of 3.7 m and stem diameters of 8.6 cm) to bigger ones (heights of 11.3 m and stem diameter of 28.3 cm).This implies a subsequent variability in biomass estimations from these variables, as can be seen in standard deviation values of aboveground biomass (AGB), Wr and TB.Table 2 shows the Spearman correlation coefficients obtained for each of the ALS-derived metrics.TB is strongly correlated to the upper ALS height percentiles, particularly with the P 25 to P 60 (0.95 in all of them), being positive in all cases.The greater the height of the trees within the plot the greater the TB value as the stems and branches are more developed.The same interpretation presents the canopy height metrics, such maximum (0.88), mean (0.95) and mode elevation (0.92).In the case of metrics related to the variability of canopy height, the skewness of the height distribution of the returns, which describes the degree of concentration of heights around high or low values, i.e.,  Figure 2 shows the TB and carbon content mapping after implementation of the prediction model in a GIS environment.
The results of this study indicate that LiDAR-PNOA data can be used to estimate the TB in monospecific Pinus halepensis Mill.forest stands.The good fit of the model obtained reveals, not only the suitability of the methodology, but also the correlation between the statistics obtained from ALS point cloud and the h and dbh measured in the field.Although the RMSE above 11,000 kg/ha may seem high, it should be noted that the inventory data covered a wide range of values as can be observed in Table 1.It is important to notice that both, the R 2 and RMSE value match those obtained by other authors, such as González-Ferreiro et al. [11].In contrast to other studies, such as García et al. [25], it was not necessary to perform a logarithmic transformation of the original variables in order to achieve a better fit of the model or to meet the assumptions of the linear regression model.
ALS metrics included in the model behave logically as can be observed in Figure 1.In this regard, it is common that prediction models include ALS-derived metrics directly related to the canopy height, variables characterizing the variability, dispersion and shape of the height distribution of returns, as is the case of kurtosis, and variables related to the canopy density and tree cover, as the percentage of all returns above 1 m.
On the other hand, the no contemporaneity of the ALS data acquisition and the field campaign was considered irrelevant, as no significant changes in forest structure occurred in that time lapse.Furthermore, the low point density of the data did not affect the existence of high and significant correlations between TB derived from field data and ALS-derived metrics.
In the future, it would be advisable to analyze the possibility of estimating forest biomass disaggregated in different fractions (stem, branches, leaves and roots) in order to model the carbon storage in detail and improve forest stand management.Besides, it would be desirable to study the usefulness of the PNOA-LiDAR data to characterize other forest species in large areas, as it could increase the accuracy of the results and reduce the costs with respect to traditional inventories.

Conclusions
Active remote sensing with ALS technology brings a new perspective to forest inventories providing 3-D information of the vegetation structure, as well as an improvement in the accuracy of the results and a reduction in costs.In this regard, the benefit of estimating the TB is even bigger, as with this methodology systematic sampling destructive procedures are not necessary.This work has proved the usefulness to estimate TB and carbon content in monospecific Aleppo pine stands representative of the continental Mediterranean climate.In the light of the results obtained, the adopted methodology was suitable.Our results match those achieved in other types of forests.Considering that PNOA project will presumably produce more LiDAR products in the next years, it would be desirable to evaluate the adequacy of the generated model to future captures of data and to develop new models of biomass components as well as in other types of Mediterranean forests.

Figure 1 .
Figure 1.Metrics associated to the vertical distribution of ALS returns in three selected field plots representative of the study area: (a) Low height, (b) medium height and (c) great height Aleppo pines.

Figure 2 .
Figure 2. Map of the TB and carbon content in two selected areas of the study site.(a) Heterogeneous Aleppo pine forest with lower stand density.(b) Includes two distinct homogeneous areas, differentiated by stand density and canopy height.

Table 1 .
Summary of the field plot data (n=45).

Table 3 .
Model summary for estimating TB.