Prediction of Aboveground Biomass from Low-Density LiDAR Data: Validation over P. radiata Data from a Region North of Spain

Estimation of forestry aboveground biomass (AGB) by means of aerial Light Detection and Ranging (LiDAR) data uses high-density point sampling data obtained in dedicated flights, which are often too costly for available research budgets. In this paper we exploit already existing public low-density LiDAR data obtained for other purposes, such as cartography. The challenge is to show that such low-density data allows accurate biomass estimation. We demonstrate the approach on data available from plantations of Pinus radiata in the Arratia-Nervión region, located in Biscay province located in the North of Spain. We use public data gathered from the low-density (0.5 pulse/m2) LiDAR flight conducted by the Basque Government in 2012 for cartographic production. We propose a linear regression model based on explanatory variables obtained from the LiDAR point cloud data. We calibrate the model using field data from the Fourth National Forest Inventory (NFI4), including the selection of the optimal model variables. The results revealed that the best model depends on two variables extracted from LiDAR data: One directly related with tree height and a second parameter with the canopy density. The model explained 80% of its variability with a standard error of 0.25 ton/ha in logarithmic units. We validate the predictions against the biomass measurements provided by the government institutions, obtaining a difference of 8%. The proposed approach would allow the exploitation of the periodic available low-density LiDAR data, collected with territorial and cartographic purposes, for a more frequent and less expensive control of the forestry biomass.


Introduction
Assessing forest resources has gained increased attention by governments worldwide in the last few decades due to increased awareness about global climate change, and a greater appreciation of the ecosystem services provided by forests [1]. Forests play a dual role in the global carbon cycle: (i) As an important carbon sink removing carbon dioxide through photosynthesis and converting the photosynthate to forest biomass; and (ii) as a carbon source by releasing carbon dioxide through respiration, wildfires, and decomposition [2]. Thus, there is great uncertainty over whether forests will be a carbon sink or source in the future. If forests are well managed, and timber is used for long-term products such as buildings, then forest management will result in a net reduction of atmospheric (average height, dominant height, or mean diameter) [13][14][15], even allowing the discrimination between tree species [16]. Forest biomass can be estimated on the basis of these variables by diverse approaches [17][18][19][20], which can be categorized depending on the footprint size and the object under study (plot size or individual tree). Early studies used small-footprint (discrete) LiDAR data to estimate biophysical properties of forest stands at plot size. They started in the 1980s using profiling lasers [21]. Later approaches measure the strength of the return of the laser pulse at each sample point, which depends on the surface reflectivity. The increase of the pulse frequency in modern LiDAR systems provides an increase of pulse density, allowing to detect treetops, gaps between the forest crowns, and individual tree analysis [22,23]. Full-waveform systems capturing the returned energy in an almost continuous fashion have been used highly accurate biomass estimation [24,25].
Specific LiDAR data capture campaigns for biomass measurement are very expensive, hindering the general application of the technology in forestry management. However, this obstacle can be overcome. Some institutions carry out periodic LiDAR data capture campaigns to build digital terrain and surface models, mostly for cartographic purposes. These datasets are published by the administration institutions, such as the Basque Government, and can be exploited for other applications [26]. The main disadvantage of these datasets are their low sampling density (0.5 pulse/m 2 ).
The purpose of this paper is to show that these low-density datasets can be useful for forest biomass estimation; therefore, allowing for their systematic use in forest management. We develop a linear regression model to estimate aboveground biomass (AGB) after optimal selection of LiDAR data features using the HAZI's allometric equations for P. radiata applied over the NFI4 data for validation purposes. The specific dataset used to demonstrate our approach is a low sampling density LiDAR dataset covering the whole area of Spain collected in the National Plan of Aerial Orthophotography (PNOA) carried out in 2011.

Study Area
The study area is the Arratia-Nervión region, located in the province of Biscay (Basque Country), in the northern part of Spain ( Figure 1). This region encompasses 14 municipalities, covering a total area of 400 km 2 . The average altitude of the region is 465 m, with an average slope of 18.6°. High slope (30-45°) areas are frequent across the entire region. Pine forests of P. radiata are the most important land cover in the Basque Country, 125,000 ha, accounting for 32% of the forested area in the Basque Country, equivalent to 49% of the area covered by this species in Spain. These pine forests provide 23% of the volume of large trees in the Basque Country, and 44% of the volume including bark. Forests of P. radiata are cultivated at altitudes below 600 m.
According to the NFI4, in Arratia-Nervión 16,260 ha out of 28,065 ha of forest areas belong to the P. radiata D. Don. tree species, representing over 60% of the tree specimens of the area (Figure 2). Coniferous plantations have progressively replaced the native tree species. The only invading species is the Robinia pseudoacacia with a minor presence of 17 ha in the entire region. Native species, such as Quercus ilex or Fagus sylvatica, represent only 12.5% (3511 ha) of the total forested areas in the region.

Field Data Collection
The dendrometric data used in this study were collected for the NFI4 in the Basque Country between 17 January and 15 June 2011. The circular plots were placed at the vertices of the UTM cartographic system (European Datum 1950) kilometric grid in areas classified as forest, locating a plot for every square kilometer ( Figure 3). The trees in the plots were measured using the methodology defined by the Nature Conservation Institute [27] based on nested plots, where each plot is subdivided into four circular plots of variable radius 5, 10, 15, and 25 m, representing a maximum area of approximately 0.2 ha for the biggest radius. The nested plot method is suitable when variability in the tree diameter exists [28]. The trees were measured depending on their Diameter at Breast Height (DBH). Minimum diameter for inclusion in the measurement was 7.5 cm. When tree DBH was between 7.5 and 12.5 cm it was included in the 5 m radius subplot, between 12.5 and 22.5 cm in the 10 m radius subplot; and between 22 and 42.5 cm in the 15 m radius subplot.

Field Data Collection
The dendrometric data used in this study were collected for the NFI4 in the Basque Country between 17 January and 15 June 2011. The circular plots were placed at the vertices of the UTM cartographic system (European Datum 1950) kilometric grid in areas classified as forest, locating a plot for every square kilometer ( Figure 3). The trees in the plots were measured using the methodology defined by the Nature Conservation Institute [27] based on nested plots, where each plot is subdivided into four circular plots of variable radius 5, 10, 15, and 25 m, representing a maximum area of approximately 0.2 ha for the biggest radius. The nested plot method is suitable when variability in the tree diameter exists [28]. The trees were measured depending on their Diameter at Breast Height (DBH). Minimum diameter for inclusion in the measurement was 7.5 cm. When tree DBH was between 7.5 and 12.5 cm it was included in the 5 m radius subplot, between 12.5 and 22.5 cm in the 10 m radius subplot; and between 22 and 42.5 cm in the 15 m radius subplot. Finally, the largest trees with DBH higher than 42.5 cm were included in a 25 m radius plot ( Figure 3). To obtain the values per ha depending on the plot radius, an expansion factor was applied (Equation (1)): where f is the expansion factor and R is the plot radius. Finally, the largest trees with DBH higher than 42.5 cm were included in a 25 m radius plot ( Figure  3). To obtain the values per ha depending on the plot radius, an expansion factor was applied (Equation (1)): ( where f is the expansion factor and R is the plot radius. We selected plots from the 118 plots in the study region where the dominant species was P. radiata (i.e., the percentage of P. radiata specimens was above 80%). Point clouds corresponding to the selected plots were compared with their corresponding orthophotos in order to detect differences that could interfere in the final result. We discarded plots where silvicultural treatments had decreased the population of P. radiata significantly, and plots with obvious mistakes in the point cloud classification. After this data cleansing, 55 plots remained ( Figure 4).  We selected plots from the 118 plots in the study region where the dominant species was P. radiata (i.e., the percentage of P. radiata specimens was above 80%). Point clouds corresponding to the selected plots were compared with their corresponding orthophotos in order to detect differences that could interfere in the final result. We discarded plots where silvicultural treatments had decreased the population of P. radiata significantly, and plots with obvious mistakes in the point cloud classification. After this data cleansing, 55 plots remained ( Figure 4). The tree diameter at breast height (DBH mm) of the tree (1.3 m) was measured using a caliper in two perpendicular directions. The tree height (m) was determined by a hypsometer. These two measures are the independent variables of the allometric equations. In our data sample, the minimum, mean, and maximum values for the tree diameter and height are 10.  The tree diameter at breast height (DBH mm) of the tree (1.3 m) was measured using a caliper in two perpendicular directions. The tree height (m) was determined by a hypsometer. These two measures are the independent variables of the allometric equations. In our data sample, the minimum, mean, and maximum values for the tree diameter and height are 10.   The tree diameter at breast height (DBH mm) of the tree (1.3 m) was measured using a caliper in two perpendicular directions. The tree height (m) was determined by a hypsometer. These two measures are the independent variables of the allometric equations. In our data sample, the minimum, mean, and maximum values for the tree diameter and height are 10.

Field Plot Positioning
The positioning of the plots was carried out, according to the NFI4 specifications, with a GPS navigator obtaining autonomous observations without any differential correction. In order to study the effect of the plot positioning error, we have shifted the plot position by 10 m in each of the compass rose directions, obtaining eight new plots around the original one as shown in Figure 6. We assume that the plots are embedded in a homogeneous forest area, so that these shifts will have no effect.

Field Plot Positioning
The positioning of the plots was carried out, according to the NFI4 specifications, with a GPS navigator obtaining autonomous observations without any differential correction. In order to study the effect of the plot positioning error, we have shifted the plot position by 10 m in each of the compass rose directions, obtaining eight new plots around the original one as shown in Figure 6. We assume that the plots are embedded in a homogeneous forest area, so that these shifts will have no effect.
We apply the Cohen's kappa concordance test measuring the agreement between two observers [29] to evaluate the similarity between the nine samples. In this case, we have nine "observers" corresponding to the 95% percentile of the height obtained from the shifted nine samples. These agreement values range from 0 to 1, where 1 implies that the compared metrics are identical, and 0 that there is no agreement. Figure 6. Representation of the simulated translation of the plot 443. Each color represents the gathered LiDAR data for the mentioned nine plots (ETRS89 UTM zone 30 North reference system).

LiDAR Data
The LiDAR data were acquired during the summer of 2012 between 12 July and 28 August. The entire Basque Autonomous Community area was flown over using a Lite Mapper 6800 Airborne Laser Scanner at an average altitude aboveground level and average speed of 1100 m and 67 m/s, respectively. The pulse repetition frequency was 100 kHz and the scan frequency was 70 kHz. The maximum scan angle was 60° with a beam divergence <0.5 mrad. The average point density was 0.5 points/m 2 . The spatial localization of the points was obtained with a precision <10 cm. The reference system is the European Terrestrial Reference System 89 (ETRS89) and the coordinate system is UTM Figure 6. Representation of the simulated translation of the plot 443. Each color represents the gathered LiDAR data for the mentioned nine plots (ETRS89 UTM zone 30 North reference system).
We apply the Cohen's kappa concordance test measuring the agreement between two observers [29] to evaluate the similarity between the nine samples. In this case, we have nine "observers" corresponding to the 95% percentile of the height obtained from the shifted nine samples. These agreement values range from 0 to 1, where 1 implies that the compared metrics are identical, and 0 that there is no agreement.

LiDAR Data
The LiDAR data were acquired during the summer of 2012 between 12 July and 28 August. The entire Basque Autonomous Community area was flown over using a Lite Mapper 6800 Airborne Laser Scanner at an average altitude aboveground level and average speed of 1100 m and 67 m/s, respectively. The pulse repetition frequency was 100 kHz and the scan frequency was 70 kHz. The maximum scan angle was 60 • with a beam divergence <0.5 mrad. The average point density was 0.5 points/m 2 . The spatial localization of the points was obtained with a precision <10 cm. The reference system is the European Terrestrial Reference System 89 (ETRS89) and the coordinate system is UTM for the thirtieth time zone. The dataset was divided into sheets of 2 × 2 km of extension, classified into eight classes: Unclassified, Ground, Low Vegetation, Medium Vegetation, High Vegetation, Building, Low Points, and Reserved. The data are publicly available at: ftp: //ftp.geo.euskadi.eus/lidar/LIDAR_2012_ETRS89/LAS/.

Orthophotos
The flight campaign carried out by the Basque Government from 23 July to 28 August 2012 produced orthophotos with a spatial resolution of 25 cm/pixel, which were used to detect possible defects in the NFI4 data, and contradictions between NFI4 and LiDAR data. These orthophotos were downloaded from the Spatial Data Infrastructure (SDI) of the Basque Country Government from the following site: ftp://ftp.geo.euskadi.eus/cartografia/Cartografia_Basica/Ortofotos/ORTO_2012/Hojas_JPG/5000/.

Biomass Field Calculation
Volume per tree was calculated using an allometric model developed by the HAZI Foundation [30] based on the destructive sample of 732 P. radiata specimens extracted from locations distributed across the Basque region. The sample was carried out between the years 1990 and 2001. HAZI's model uses the diameter at breast height (d mm) and total tree height (h m) as input variables (Equation (2)): The resulting biomass value was corrected adding 4.04% of the obtained volume, in order to account for the tree branches and upper part of the trunk discarded for wood production reasons. This fixed correction quantity is specified per species, independent of the characteristics of the forest stand. The biomass of each sample plot was computed applying the allometric equation (Equation (2)) to the measured trees in the plot, adding them and computing the extrapolation to the size of control plot (25 m radius) using the expansion factor defined in Equation (1). Stand volume is obtained as the addition of the expanded values of the tree volumes of each plot.

LiDAR Data Processing and Overall Process
The complete LiDAR processing workflow is displayed in Figure 7a. In the initial step, the original LiDAR data (stored in LAS format files) were cleaned to identify and filter out suspected erroneous points. This is an important step, as these points can introduce errors in ensuing calculations. For this task, the altimetry value range of the point cloud was divided into equal size (15 m) intervals, counting the number of returns falling in each interval. Automatic detection of erroneous returns is based on the detection of empty layers between non-empty layers. We filtered out 0.06% of the total number of echoes. In the following steps we use the FUSION/LDV (LiDAR Data Viewer) (http://forsys.cfr.washington.edu/fusion/fusionlatest.html) free software. The next step was the selection of the points corresponding to the circular plots of 25 m radius considered in the NFI4, maintaining the dimensions of the control plots.
The returns classified as ground were used to create a digital terrain model (DTM) with a spatial resolution of 1 m 2 /pixel. A digital surface model (DSM) was created with the highest returns of the point cloud over the sample plots, assigning the elevation of the highest return within each grid cell to the grid cell center. A canopy height model (CHM) was obtained by subtracting the DTM from the DSM. The CHM characterizes the tree canopy and it is able to give the canopy height directly. For each plot a large collection of metrics, that can be used as regressors for biomass prediction, was computed over all returns above a 2 m threshold [31,32], including height distributions, canopy density metrics, and descriptive statistics (these metrics are enumerated in Appendix A). The returns classified as ground were used to create a digital terrain model (DTM) with a spatial resolution of 1 m 2 /pixel. A digital surface model (DSM) was created with the highest returns of the point cloud over the sample plots, assigning the elevation of the highest return within each grid cell to the grid cell center. A canopy height model (CHM) was obtained by subtracting the DTM from the DSM. The CHM characterizes the tree canopy and it is able to give the canopy height directly. For  Canopy density metrics have demonstrated their usefulness as predictors of forest parameters such as mean height, dominant height, mean diameter, basal area, or timber volume [15]. We used the PostGis environment to obtain new metrics related to the canopy point density for each sample plot. The point cloud and the DTM were transformed into tables using routines developed with the PostGis spatial database extender for the PostgreSQL Database Management System. We divided the point cloud into 10 vertical layers of equal height. We set the lower limit at 2 m from the ground, to avoid shrubs, and the upper limit as the 95% percentile of the distribution of heights. This last metric was chosen instead of the maximum height, due to the stability demonstrated in previous studies [33,34]. Then, the routine counts the fraction of points falling inside each layer. That way, 10 canopy densities were computed (denoted as tr_1,..., tr_10). Figure 7b provides a diagram describing the overall process carried out in the study. From the LiDAR data, we extract the LiDAR features that are the regressors of the regression model, we carry out a feature selection on the basis of the predictive performance and posthoc statistical tests of the regression model, selecting an optimal subset of LiDAR features. The NFI4 field data is used as the ground truth for model training and validation, which includes a cross-validation experiment as well as direct comparison of the predicted biomass with the NFI4 data for a site not included in the training data. Finally, we make additional comparison with the estimations provided by HAZI institute, which are produced by an allometric equation whose input variables are extracted from the LiDAR data.

Regression Analysis
This method computes a prediction of the variable under study as a linear combination of a set of regressor variables, often called features or input factors (Equation (3)): Although it can be argued that the constant offset b0 does introduce a bias in the model, because setting all coefficients to zero would lead to a non-zero result, we are pretty sure that the specific context of the modeling in this paper does not include such degenerate case. We have extracted a large number of metrics from LiDAR data which can be used as regressors for biomass prediction. Some of them can be redundant or irrelevant. In order to select the most informative metrics to be used as regressors, firstly the Pearson correlation coefficient (r) was calculated between the biomass values and our LiDAR metrics. Secondly, the same operation was carried out between the logarithm of biomass values and our LiDAR metrics in order to check the linear relation. Dispersion graphics were used to ascertain the linear relation between the biomass and each of the candidate regressors.
Subsequently, we explored the quality of the prediction using all possible combinations of feature selections over the LiDAR metrics up to three variables per model. The adjusted coefficient of determination (R 2 adj) was used to assess the quality of the adjustment. R 2 adj represents the proportion of the variability explained by the adjusted model [35] after application of the correction factor described below. Additionally, standard error of the estimations (SEEs) were calculated in order to compare with other studies results. Akaike information criterion (AIC) is an index based on in-sample fit that can be used as an estimate of the likelihood of a model to predict the future values [36]. Our optimal feature selection corresponds to the model that has minimum AIC.
To check the hypothesis of the linear regression technique, several tests were applied to the trained models: Shapiro-Wilk test to verify the normality of the residuals, Breusch-Pagan test to analyze the homoscedasticity of the residuals, Durbin-Watson test to detect possible dependencies between the data, Variance Inflation Factor (VIF) to detect collinearity problems in the model [37], Ramsey´s RESET linearity test to verify lineal relations, and, finally, Bonferroni´s test to find statistically significant atypical values. All tests, except VIF, were calculated at the 95% level of significance.
When using the log model, it is necessary to compute the inverse of the logarithmic transformation, which may introduce some bias in the distribution of the estimated biomass values leading to under-estimation of the biomass [38]. To minimize the bias introduced in the model, a correction factor was applied, which depended on the standard error of the estimate (SEE). Once the SEE is calculated, the correction factor is computed as follows: We carried out a sensitivity analysis of the model using an extended FAST (Fourier amplitude sensitivity test), because it is able to capture the influence of each regressor over its full range of variation. The objective of the sensitivity analysis (SA) of a regression model output is to ascertain how its output depends on its regressors [39]. This test allows the computation of the total contribution of each regressor to the variance of the output, as well as all the contribution of the interaction terms involving that regressor.
Let the fitted regression model be denoted as: where X i is a regressor model. The, S i is the contribution of regressor X i to the variance of the output Y as specified in the following: where E(Y|X i ) is the expectation of the output Y conditioned to input factor X i and V(Y) the unconditional variance of the model output, which can be decomposed into the conditional regressor variances V i.. as follows: Dividing both sides of the equation above by V(Y) we obtain: where S i is the first-order sensitivity index for regressor X i , and S ij is the second-order sensitivity index for the interaction between regressors X i and X j with j i. The total contribution of regressor X i to output Y variability is computed as follows: S Ti is the total effect on the output variation due to factor X i , adding its first-order effect and all higher-order effects due to interactions with other regressors. When the sum of the first-order index and total-effect index of a variable is not equal to one, interactions among factors in the model may occur. Additionally, we compute the fraction of the output variance arising from the uncertainty of each regressor X i , and the complementary set of regressors, denoted D 1 and D t , respectively.
The biomass and the LiDAR data of one of the municipalities included in the study region, Orozko (Figure 8), was used to fit the LiDAR regression model. This municipality encompasses almost the 24% of the total population of P. radiata in the study area. More than 4260 ha of P. radiata were used for the model parameter estimation, taking into account all the polygons of the species located in the study area according to the NFI4 for the Basque Country and their occupation percentage. The biomass estimations provided by HAZI foundation for this municipality were used as the ground truth biomass. For its estimation, HAZI used NFI4 data and LiDAR data from the flight of 2012. They build a linear regression model of the wood volume using the mean height of the LiDAR points above 4 m aboveground as the single regressor of the model. Then, the biomass was calculated by adding, to the volume, 4% of the obtained volume, because when calculating the volume, the tree's branches and a part of the trunk is not taken into account because of wood production reasons. This fixed quantity is stipulated per species, independent of the characteristics of the forest stand.
Forests 2018, 9, x FOR PEER REVIEW 12 of 26 of 2012. They build a linear regression model of the wood volume using the mean height of the LiDAR points above 4 m aboveground as the single regressor of the model. Then, the biomass was calculated by adding, to the volume, 4% of the obtained volume, because when calculating the volume, the tree's branches and a part of the trunk is not taken into account because of wood production reasons. This fixed quantity is stipulated per species, independent of the characteristics of the forest stand.

Validation
The validation process consists in the comparison between the biomass predictions and the observations over a set of measures that are different to the ones used in the model adjustment [40]. For this purpose, we used the data from Encartaciones region (Biscay) because of its large amount of P. radiata forests. This region is composed of eighteen municipalities, with a total area of 55,300 ha. Specifically, we use the data from the Gordexola municipality to validate the model because it includes almost 20% of the entire population of P. radiata of the region (Figure 15). More than 2600 ha of P. radiata were used to validate the model, according to the NFI4 polygons ( Figure 8).
Besides, a k-fold cross-validation technique was applied with the same data that was used for the regressor selection, whereby the entire number of sample plots was divided into k = 10 subsamples composed of similar number of plots. The parameter estimation and validation was repeated k times, taking, each time, a different fold as the test and the remaining k-1 folds as the data for model training. The mean square error (MSE), was calculated as follows (Equation (10)): where, lnYi is the natural logarithm of the values of the dependent variable, l n i Y is the natural logarithm of the model estimations, N the number of cases of the sample. In fact, we will use the root of the MSE (RMSE) [41].

Validation
The validation process consists in the comparison between the biomass predictions and the observations over a set of measures that are different to the ones used in the model adjustment [40]. For this purpose, we used the data from Encartaciones region (Biscay) because of its large amount of P. radiata forests. This region is composed of eighteen municipalities, with a total area of 55,300 ha. Specifically, we use the data from the Gordexola municipality to validate the model because it includes almost 20% of the entire population of P. radiata of the region (Figure 15). More than 2600 ha of P. radiata were used to validate the model, according to the NFI4 polygons ( Figure 8).
Besides, a k-fold cross-validation technique was applied with the same data that was used for the regressor selection, whereby the entire number of sample plots was divided into k = 10 subsamples composed of similar number of plots. The parameter estimation and validation was repeated k times, taking, each time, a different fold as the test and the remaining k-1 folds as the data for model training. The mean square error (MSE), was calculated as follows (Equation (10)): where, lnY i is the natural logarithm of the values of the dependent variable, lnŶ i is the natural logarithm of the model estimations, N the number of cases of the sample. In fact, we will use the root of the MSE (RMSE) [41].

Results
Firstly, we study positioning error effect, reporting the results of the Cohen concordance test. Table 1 shows the agreement between the original plot and the eight translated plots for plot 443 (the same procedure was applied to the remaining 54 plots). The agreement level among plots shown in the upper triangle of the matrix is close to 1, with a minimum value of 0.9. Hence, the eight plots created around the original one are similar enough to assure the low influence of the positioning error.
We computed the correlation matrix between the LiDAR data based regressors (specified in detail in Table A1) and the field biomass values and their logarithm. The highest correlated regressors were the height percentiles reaching maximal correlations of r = 0.80 and r = 0.86 with the biomass and the log-biomass, respectively. Although the regressors induced from canopy density measures were not highly correlated with the biomass, their relation to the biomass was statistically significant for all of them, reaching a maximum r = 0.62.
The linear relation can be assessed using dispersion diagrams in Figure 9 [42]. The shaded band is a pointwise 95% confidence interval on the fitted values (the blue line).

Results
Firstly, we study positioning error effect, reporting the results of the Cohen concordance test. Table  1 shows the agreement between the original plot and the eight translated plots for plot 443 (the same procedure was applied to the remaining 54 plots). The agreement level among plots shown in the upper triangle of the matrix is close to 1, with a minimum value of 0.9. Hence, the eight plots created around the original one are similar enough to assure the low influence of the positioning error.
We computed the correlation matrix between the LiDAR data based regressors (specified in detail in Table A1) and the field biomass values and their logarithm. The highest correlated regressors were the height percentiles reaching maximal correlations of r = 0.80 and r = 0.86 with the biomass and the log-biomass, respectively. Although the regressors induced from canopy density measures were not highly correlated with the biomass, their relation to the biomass was statistically significant for all of them, reaching a maximum r = 0.62.
The linear relation can be assessed using dispersion diagrams in Figure 9 [42]. The shaded band is a pointwise 95% confidence interval on the fitted values (the blue line).   After the identification of the highest log-biomass correlated variables, and assessing their relationship, we fitted the best models using only one variable, two variables, or the combination of three variables.
The inclusion of three variables improved the adjusted coefficient of determination (R 2 adj = 0.81) in some cases, but introducing the third variable was not statistically significant in the models, hence, we discarded using more than two variables. Since the models presented in the table have very similar R 2 values, further statistical analysis was undertaken to decide which model better fulfilled the assumptions of the linear regression analysis.
As can be seen in Table 2, the ten best fitting MLR models obtained very similar results in the inference tests. The models show identical numerical values for the R 2 adj value (0.79) and standard error (0.25 ton/ha in logarithmic units), the remaining columns reporting tests results had similar values too, including the detection of outliers according to Bonferroni´s test in all the models, where no outliers were detected. Regarding the normality of the residuals, it is not possible to reject this hypothesis in any model. The results are slightly better in the models using the 99% percentile of the height. In the case of the Breusch-Pagan test, the p-values values do not vary too much among regressor subsets, so it will be acceptable for all the models in the table, corresponding the best results (in the sense of non-rejection of null hypothesis) to the eighth model. The values obtained for the Durbin-Watson test and the variance inflation factor are very similar for all the models, concluding that no autocorrelation or collinearity problems were detected. The RESET linearity test results confirm the null hypothesis of a linear functional dependence of the biomass on the regressors for all tested models. Table 2. Accuracy values and test p-values obtained for the ten best models; p-values are shown for all the tests except VIF. p99, p95, and p90 are the 99%, 95%, and 90% percentiles of the laser canopy heights, respectively; tr_2, tr_3, and tr_4 are the canopy densities corresponding to the second, third, and fourth layers, respectively; allabovemean = (all returns above mean height)/(total returns). R 2 adj = adjusted coefficient of determination; SE = Standard Error; AIC = Akaike information criterion; SW = Shapiro-Wilk residuals normality test; BP = Breustch-Pagan residuals homoscedasticity test; DW = residuals autocorrelationi test; VIF = Variance Inflation Factor; reset = Ramsey´s RESET linearity test; B = Bonferroni outlier test. For the optimal selection of variables, we focus on the Breusch-Pagan test as heteroscedasticity is to be avoided in regression models, thus the selected model is {p95, tr_3} (i.e., the one composed of the 95% percentile of the LiDAR measurement of canopy heights (p95), and the percentage of points above the third layer from the total number of returns (tr_3)). The equation corresponding to this model is the following (Equation (11)):
As it has been commented in Section 2.6.3, the logarithmic values of the estimated biomass must be transformed to arithmetic values. In this case the value obtained for the correction factor was 1.032, hence, the biomass equation reads: .Biomass (ton/ha) = 1.032 · e (3.77418+(0.06729·p95)+(0.54792·tr_3)) (12) Extensive diagnostic graphics, shown in Figures 10-12, were carried out on the best model. It does not show heteroscedasticity behavior, as it can be checked upon inspection of the "residual vs. fitted" values plot, and the "Scale Location" plot, where no tendency of the residuals can be identified. The Quantile-Quantile (Q-Q) plot confirms the normality assumption of the residuals, previously checked by a Shapiro-Wilks test. There are "light tails" in the Q-Q plot, more pronounced in the right side of the plot, so it is detecting that there were more extreme values than would be expected for a truly normal distribution. The results of the Bonferroni test already removed doubts over the existence of outliers or atypical values. The dispersion plots in Figures 11 and 12 confirm the linear dependency hypothesis, and the normality of the residuals of the model predictions.  (12) Extensive diagnostic graphics, shown in Figures 10-12, were carried out on the best model. It does not show heteroscedasticity behavior, as it can be checked upon inspection of the "residual vs. fitted" values plot, and the "Scale Location" plot, where no tendency of the residuals can be identified. The Quantile-Quantile (Q-Q) plot confirms the normality assumption of the residuals, previously checked by a Shapiro-Wilks test. There are "light tails" in the Q-Q plot, more pronounced in the right side of the plot, so it is detecting that there were more extreme values than would be expected for a truly normal distribution. The results of the Bonferroni test already removed doubts over the existence of outliers or atypical values. The dispersion plots in Figures 11 and 12 confirm the linear dependency hypothesis, and the normality of the residuals of the model predictions.    The minimum, mean, and the maximum estimated biomass per plot was 69.06, 294.77, and 611.21 ton/ha, respectively. As can be shown in Figure 13, the distribution of the biomass estimations does not appear to be normal. We visualize in Figure 14 the spatial location of the sample plots and the magnitude of the biomass estimated at each plot.   The minimum, mean, and the maximum estimated biomass per plot was 69.06, 294.77, and 611.21 ton/ha, respectively. As can be shown in Figure 13, the distribution of the biomass estimations does not appear to be normal. We visualize in Figure 14 the spatial location of the sample plots and the magnitude of the biomass estimated at each plot. The minimum, mean, and the maximum estimated biomass per plot was 69.06, 294.77, and 611.21 ton/ha, respectively. As can be shown in Figure 13, the distribution of the biomass estimations does not appear to be normal. We visualize in Figure 14 the spatial location of the sample plots and the magnitude of the biomass estimated at each plot.  We carry out statistical inference over the linear regression model with optimal regressor selection obtaining the upper and lower bounds shown in Table 3. The average biomass estimate over   We carry out statistical inference over the linear regression model with optimal regressor selection obtaining the upper and lower bounds shown in Table 3. The average biomass estimate over We carry out statistical inference over the linear regression model with optimal regressor selection obtaining the upper and lower bounds shown in Table 3. The average biomass estimate over the 55 sample plots is 294.77 ton/ha. Its lower and upper bounds of the 95% confidence interval are 189.79 and 460.37 ton/ha, respectively. The average inference over the biomass estimated in the 55 plots is 294.77 ton/ha, establishing as low and upper limit at a 95% confidence interval, 189.79 and 459.05 ton/ha, respectively.
In addition to these values, simulations in each of the 55 plots were performed, doing 10,000,000 iterations in each plot, at a 95% confidence interval applying Monte Carlo simulation. Monte Carlo uses random samples to simulate data for a certain mathematical model. The results were very similar to the ones previously obtained, in this case, the average value of the biomass is exactly the same and the low and upper limit are 190.35 and 460.37 ton/ha respectively.
Sensitivity analysis with the extended FAST test [43] was carried out, using M = 4 as interference factor and n = 17 as sample size. The results of the analysis are displayed in Table 4. The model appeared to be non-additive because the sum of the first-order index and the total-effect index is not equal to one for both regressors. Hence, there must be interactions among them. The existence of interactions may imply, for instance, that extreme values of the output Y are uniquely associated with particular combinations of model inputs, in a way that is not described by the first-order effects (S i ). The total-effect indices (ST i ) show that only p95 is taking part in the interactions, because ST i > S i . D 1 represents the portion of the output variance arising from the uncertainty of factor i, while D T is the variance for the complementary set D(-i) [44]. First-order effect (S i ) is bigger for the density metric than for the percentile of the height, which implies that this variable is more sensible in the model. D 1 and D T corroborates this tendency, pointing that the density metric contributes more than the height percentile to the output variance, because D 1 is bigger for the density metric than for the height percentile. Table 4. Sensitivity analysis results. Si = first-order effects; STi = total-effect index; D 1 = the portion of the output variance arising from the uncertainty of factor i; D T = variance for the complementary set D(-i).

Application of the Selected Model
The fitted optimal regression model was applied to LiDAR data extracted from the Orozko municipality, in Arratia-Nervión, comparing its biomass estimation with the biomass estimation given by the HAZI foundation using allometric equations applied to field measurements, finding a difference of 8%, 992,747.53 vs. 915,851.29 ton. The spatial distribution of the estimation of our regression model is visualized in Figure 15.

Validation of the Selected Model
As detailed in Section 2.6.4, the first validation step applied the developed model in the Gordexola municipality (an independent area). A difference of 121,734.65 ton was obtained between the biomass given by the HAZI foundation (561,868.93 ton) and the one estimated in this study (683,603.58 ton), representing approximately 21% of the total biomass of the Gordexola municipality.
In a second step, a k-fold cross-validation was also applied, in which the initial sample, composed of the same 55 plots used to estimate the model, was subdivided in five equal folds and the algorithm fitted to the model using the k-1 remaining folds. Eleven elements were used per fold, with the mean value of the addition of the square of the residuals oscillating between 0.0674 and 0.0682 logarithmic units, more or less 0.7% of the mean value of the field biomass, suggesting that the prediction values made were quite accurate.

Discussion
We proposed and tested a linear regression model of the logarithm of the biomass using regressors extracted from low point density LiDAR data. A ground truth biomass estimate was constructed by the application of an allometric model with dasometric values from the NFI4. Our optimal model has only two independent variables: a height percentile and a canopy density metric, both obtained from the LiDAR data. The model explains 80% of the variance of the sample with an RMSE of 0.5 in logarithmic units. Our results presented here are within the range of those reported by similar studies using LiDAR technology to characterize AGB in coniferous forest, profiting from higher density sampling LiDAR data [32,[45][46][47]. Similar results to ours (R 2 = 0.74; SE = 0.2) were reported by Hall et al. [48] in an area in Colorado, USA, for Pinus ponderosa and Pseudotsuga menziesii species. The density points in the American study were higher than the one used in our study (1.23 points/m 2 ), but in their study the only variable introduced in the exponential model was a density metric, no metrics derived for the tree height were included.

Validation of the Selected Model
As detailed in Section 2.6.4, the first validation step applied the developed model in the Gordexola municipality (an independent area). A difference of 121,734.65 ton was obtained between the biomass given by the HAZI foundation (561,868.93 ton) and the one estimated in this study (683,603.58 ton), representing approximately 21% of the total biomass of the Gordexola municipality.
In a second step, a k-fold cross-validation was also applied, in which the initial sample, composed of the same 55 plots used to estimate the model, was subdivided in five equal folds and the algorithm fitted to the model using the k-1 remaining folds. Eleven elements were used per fold, with the mean value of the addition of the square of the residuals oscillating between 0.0674 and 0.0682 logarithmic units, more or less 0.7% of the mean value of the field biomass, suggesting that the prediction values made were quite accurate.

Discussion
We proposed and tested a linear regression model of the logarithm of the biomass using regressors extracted from low point density LiDAR data. A ground truth biomass estimate was constructed by the application of an allometric model with dasometric values from the NFI4. Our optimal model has only two independent variables: a height percentile and a canopy density metric, both obtained from the LiDAR data. The model explains 80% of the variance of the sample with an RMSE of 0.5 in logarithmic units. Our results presented here are within the range of those reported by similar studies using LiDAR technology to characterize AGB in coniferous forest, profiting from higher density sampling LiDAR data [32,[45][46][47]. Similar results to ours (R 2 = 0.74; SE = 0.2) were reported by Hall et al. [48] in an area in Colorado, USA, for Pinus ponderosa and Pseudotsuga menziesii species. The density points in the American study were higher than the one used in our study (1.23 points/m 2 ), but in their study the only variable introduced in the exponential model was a density metric, no metrics derived for the tree height were included. Naesset and Gobakken [32] obtained higher R 2 values (0.88) than us with a very similar value for RMSE (0.25 in logarithmic units) in a Norwegian area composed primarily of Pinus sylvestris and Picea abies using a point density up to 1.2 points/m 2 . Their chosen regressors were a high percentile of the LiDAR measured tree height (90% percentile) and a low-density metric. The density metric was more statistically significant in the Norwegian model than in our study: They found that the partial value of R 2 for each variable was 0.61 and 0.21, respectively, while in our study the values were 0.74 and 0.15, respectively. The mean value of the estimated biomass was similar in both studies (150 ton/ha). Hence, differences in the model determination coefficient can be due to the geographical location, altitude, species composition, or site quality. Another relevant study [49] was carried out over the Taita Hills (55,000 ha), in southeast Kenya, with a pulse density of 3.1 points/m 2 . They used a boosted regression trees (BRT) technique for identifying regressors that better explain biomass distribution. Multilinear regression was used to predict aboveground biomass using LiDAR metrics (R 2 = 0.88 and RMSE = 52.9 Mg/ha) and the mean AGB was 123 ton/ha. In this case, the point density was significantly higher than in our study. It has been shown that pulse density below 1 pulse/m 2 has a negative influence on the quality of the prediction [50].
The results obtained by González-Ferrreiro et al. [46] for P. radiata in Galicia (Spain) were very similar to ours (R 2 = 0.74, RMSE = 40.469 ton/ha). The study area was situated in the north of the peninsula too, with similar conditions to those found in the Arratia-Nervión (mean biomass = 150 ton/ha in both cases), but their point density was very high (8 points/m 2 ). In the province of Zaragoza (Spain), the work reported in [26] estimated the biomass with good results (R 2 = 0.89, RMSE = 7326.12 kg/ha), but the mean and the maximum biomass measured in the field were much lower, along with the variability of the samples.
The Q-Q plot in Figure 10 corroborates the normality of the residuals, but moderate tails can be appreciated, which are slightly heavier in the right side of the figure. This shows that the normality of the residuals decreases when their values increase, which translates into biomass underestimation in relation with the field measurements. The plots with higher errors (>40%) correspond to the same timber stage, with tree diameters higher than 20 cm. Even the residuals do not show a clear tendency, indicating an overestimation between LiDAR heights and those measured in the field (max. 10 m), when the opposite was expected due to the difficulty of the LiDAR to detect the top of the trees [42]. These differences might be related with errors during the field data acquisition, the time mismatch between both data sets, and, finally, errors in the height determination from the LiDAR data. Apart from these circumstances, the usage of the circular nested-plots methodology applied in the NFI4 and the extrapolation of the data to the 25 m plot may be another source of error. No observation has exceeded the Cook distance in the "residual vs. leverage" plot, so they might not be influential observations, but it is true that some observations had high residuals.
No homoscedasticity issues were detected in the model, when the adjusted values of biomass were compared graphically with the residuals (Figure 11a), hence, it can be concluded that the variance of the residual is independent of the regressors. Furthermore, the linear adjustment seems to be the most suitable, as the linearity test confirmed. A clear relationship exists between the diameter measured in field and the estimated biomass (Figure 11b), as the values of the biomass increase, the values of the diameter are more disperse. This condition is well known in the use of statistics models for biological phenomena where high variability can appear: Equal dimensions of the trunk does not imply equal quantity of biomass, because it depends on the growing conditions of the trees [14]. As can be seen in Figure 12, the linear relationship between the logarithm of the biomass and the 95% height percentile is more evident than the relation between the explicated variable and the canopy density metric included in the model. As the regression proves, the p-value obtained was <0.001 for the percentile, and 0.039 for the density metric regressor, with a significance level of 0.05. But no other tendencies were deduced from the plots.
The main motivation of our study is to show that LiDAR data with lower point density than the densities reported in LiDAR data capture missions dedicated to biomass measurement can be used for AGB estimation. In this regard, some studies have estimated biomass using data with low point densities obtained by subsampling LiDAR data from flights especially designed for biomass estimation studies, with the aim of analyzing the influence of point densities in the final estimations. Authors report different conclusions. One study calculated the biomass in a forested area of central Spain, using LiDAR data with a final average point density of 11.4 points/m 2 , due to the acquisition of several overlapped scans [51]. They report results from nine different point density values: 0.25, 0.5, 1, 1.5, 2, 3, 4, 5, and 6 points/m 2 , obtaining increasing values for R 2 ranging from 0.69 to 0.92 as the point density increased. Another work [52] analyzed the importance of the point density for biomass estimation in three different sites across Ontario, Canada. They concluded that even with a very low point density of 1/100, automated LiDAR scan (ALS) is a feasible option for assessing AGB in vast areas of flat, lowland peat swamp forest. Additionally it has been found [50] that the accuracy of predicted forest structure metrics decreases as the pulse density decreases, remaining relatively high until low densities (e.g., 1 pulse/m 2 ). Due to the influence of all these factors, replication effects can be detected in the final AGB estimations. Magnussen et al. [53] found that a minimum point density of 1 point/m 2 was needed to reduce the replication variance in LiDAR-derived predictions. Known factors affecting extracted ALS-based predictors of forest inventory attributes are instrument errors, positional errors, and posts-capture data-processing errors. We must take into account that these factors limit the ability to replicate results of AGB estimation. However, the major conclusion of the reviewed studies is that useful AGB estimations can be obtained from low-density LiDAR data, as done in our study.
The difference obtained between the biomass estimation data provided by the HAZI foundation and our predictions may be due to the combination of various uncontrollable sources of variability along the data extraction and computational processes. Firstly, HAZI calculated the wood volume using the allometric equation applied to tree measures extracted from LiDAR data, while we estimated the biomass directly from the LiDAR data, avoiding the accumulation of errors due to the concatenation of various mathematical processes. HAZI subsequently added, to the estimated volume, a constant percentage of the obtained volume, which is species specific. Corona et al. [54] have derived a ratio estimator of the total volume of a forest, constituting an approximately unbiased estimator of the total volume of its sampling variance and the corresponding confidence interval. This method needs a very accurate georeference and co-registration of both LiDAR measurements and plot locations. In addition, in large areas, including forest stands with different characteristics (species, silvicultural treatments, age structures, etc.), stratification will be required.
Even if the developed biomass estimation model is local, the designed automatic and straightforward process to obtain the biomass estimations can be applied in different regions and for different species. The main limitation of the application of the explained methodology is the availability of public LiDAR and National Forest Inventory (NFI4) data. The LiDAR data used in this study are publicly available, continuous in time, and are collected in cycles of four/five years, more frequently than the data gathered in the NFI4. This study found that height percentiles were the most relevant height-related regressors for biomass prediction. As other authors have found [55], canopy density metrics have statistical significance in the model and, hence, were included. There are some major sources of unexplained variance of the AGB regression model [56][57][58]). The first is mis-registration of the ground plot locations to the LiDAR point cloud. As has been shown above, position errors are not influential in the results. The second is the application of allometric equations for biomass predictions over the ground measurements. In the NFI4 not all of the trees were measured, as we have explained, only some individuals were measured and then the values were extrapolated. The third error source would be related with the plot size, because of the disagreement between LIDAR and field plot measurements, over which trees or parts of trees are inside calibration plots. The biomass estimations presented in this study may be influenced by the exclusion in the ground data of overhanging tree crowns captured by LiDAR, which may lead to underestimation of the final result [8]. Finally, classification errors of the LiDAR data into levels could interfere with the results. Carrying out information fusion with other sensors [59,60] may alleviate these problems According to the quality curves 1, 2, and 3 for P. radiata species for the Basque Country developed in 1974, the estimated average growth per year is approximately 1 m. In the adjusted regression models, our regressors are the LiDAR heights from the flight of 2012, while our ground truth is the biomass measured in 2011, thus we have only one year of growth to account for. In the ground truth volume measurements, a correction of 4% was added to account for biomass discarded by the field researchers. We believe that the one-year growth missed can be accounted for by this operation.

Conclusions
This study has shown the application of an automatic and straightforward process to develop a local model to estimate biomass of P. radiata using LiDAR data from a low point density flight (0.5 points/m 2 ). The approach could be transferred to other areas, if LiDAR and forest inventory datasets are available, and could become a powerful tool for complimenting other traditional methods (e.g., NFI4), reducing significantly the high investment of time and money required by such methodologies.
Based on the encouraging results obtained in this study, we will propose machine learning techniques to improve results in future studies. Similarly, the incorporation of data from additional sensors (e.g., hyperspectral images or synthetic aperture radar) could help improve the model-based results. For future research, exploiting open data from the European Copernicus program could be a good avenue to obtain improved predictive models comprising satellite-borne Earth observation and in-situ data. A service component that combines these information sources will provide essential information for monitoring the terrestrial environment.