Temporal Transferability of Pine Forest Attributes Modeling Using Low-Density Airborne Laser Scanning Data

: This study assesses model temporal transferability using airborne laser scanning (ALS) data acquired over two different dates. Seven forest attributes (i.e. stand density, basal area, squared mean diameter, dominant diameter, tree dominant height, timber volume, and total tree biomass) were estimated using an area-based approach in Mediterranean Aleppo pine forests. Low-density ALS data were acquired in 2011 and 2016 while 147 forest inventory plots were measured in 2013, 2014, and 2016. Single-tree growth models were used to generate concomitant ﬁeld data for 2011 and 2016. A comparison of ﬁve selection techniques and ﬁve regression methods were performed to regress ﬁeld observations against ALS metrics. The selection of the best regression models ﬁtted for each stand attribute, and separately for both 2011 and 2016, was performed following an indirect approach. Model performance and temporal transferability were analyzed by extrapolating the best ﬁtted models from 2011 to 2016 and inversely from 2016 to 2011 using the direct approach. Non-parametric support vector machine with radial kernel was the best regression method with average relative % root mean square error differences of 2.13% for 2011 models and 1.58% for 2016 ones.


Introduction
Forest ecosystems provide economic and social benefits to humankind [1,2] requiring accurate tools to monitor their dynamics over time [3]. Over the last decades, optical remote sensing techniques have been used for monitoring forest changes at regional scales with the support of field surveys (e.g., [4,5]). However, airborne laser scanning (ALS) is better adapted to characterize forest structure [6] and estimate forest inventory parameters, providing accurate information to perform forest management and planning [3]. Furthermore, costs of ALS-based inventories are comparable to those associated with traditional ground-based ones [7,8]. Despite the great potential of this technology, multi-temporal ALS data have been utilized less, as the availability of two or more surveys in the same area has been limited by acquisition costs as well as by the need of temporal-concomitant field data (e.g., [3,6,9,10]). Recently, organizations, companies, and countries have made an effort to gather multi-temporal datasets in different years (e.g., [11][12][13]) allowing the estimation of biophysical properties in forested areas over time. As a result, height growth has been estimated using the single tree or the area-based approach [14][15][16][17][18][19][20]. Biomass and carbon dynamics has also been analyzed [3,9,[21][22][23][24][25][26][27], while less studies have estimated volume [17,18,21], basal area [17], and site index [14,28]. Multi-temporal data has also been used for quantifying fire-induced changes to forest structure [29], gaps presence [20,30,31] or detecting defoliation [32]. However, to the best of our knowledge, some relevant forest inventory attributes, such as stand density, squared mean diameter, and dominant diameter, have not been estimated using multi-temporal ALS data.
The use of low-density ALS data has been successfully used for estimating forestry attributes in different forest ecosystems (e.g., [33][34][35]), being also the case in Mediterranean forests of Spain (e.g., [36][37][38]). The analysis of the influence of point density on model predictions have been analyzed by several authors (e.g., [39][40][41]), who established that point density has little or no effect on predictions as statistical metrics remain stable [42]. Furthermore, Garcia et al. [43], Singh et al. [44], and Ruiz et al. [45] pointed out that low-density datasets were a viable solution at regional scales. Furthermore, the use of multi-temporal, low-point density data has only been explored in boreal ecosystems [17,28] and in temperate forests [24,27] but not in other ecosystems, such as Mediterranean ones, which are characterized by a higher heterogeneity in terms of forest structure.
Direct and indirect approaches have been proposed to model forest attributes using multi-temporal ALS data over time [26]. The direct approach adjusts one model for one point in time and estimates the inventory attribute for another point in time [3,9]. Previously, the model to be extrapolated was generated through regression methods that related a suite of ALS-metrics with ground-truth data. This approach allowed the temporal transferability of models reducing modeling time and fieldwork, as it was not needed to revisit them when the time between the ALS surveys was not large [28]. In contrast, the indirect approach fits two different models and estimates the variables for each point in time [3,9,17,24,25,27]. Several examples of the evaluation of these two approaches can be found in the literature. For example, when estimating biomass and carbon fluxes, Zhao et al. [3] and Meyer et al. [25] achieved better results with the indirect approach while Cao et al. [9], Skowronski et al. [24], and Bollandsås et al. [26] found slightly better performance of the direct approach.
These aforementioned modeling strategies sometimes face a challenge when lacking temporallyconcomitant field data to calibrate forest stand models [3]. To this end, forest growth models can serve as useful tools to calibrate forest stand variables in a specific point in time. Thus, empirical growth models have received particular attention since the beginning due to their usefulness. Nowadays state-space stand-level models [46], distribution-based models, and both individual-tree models and complex process-based eco-physiological models [47] have dramatically increased flexibility and realism to forest simulations. In this sense, individual-tree growth models are powerful tools to update stand variables to the Light Detection and Ranging (LiDAR) mission date. For example, the use of diameter at breast height (dbh) and the height growth values from general yield tables of the Spanish National Forest Inventory have been applied for estimating total tree biomass in Aleppo pine forests [36]. However, specific single-tree growth models calibrated with tree rings are more accurate, particularly for improving model consistency when working at regional scale.
Thus, the aim of this study is to assess temporal transferability of several forest attributes models by comparing direct and indirect approaches using low-density ALS datasets collected in 2011 and 2016. Seven forestry attributes (i.e. stand density, basal area, squared mean diameter, dominant diameter, tree dominant height, timber volume, and total tree biomass) are estimated at regional scale in the Mediterranean Aleppo pine forest. First, an indirect approach fits two different models for 2011 and 2016 and estimates stand attributes for each point in time using different ALS-metrics and model parameters. Secondly, a direct approach extrapolates the models fitted previously, using the same variables and model parameters, to the other points in time.
Furthermore, the following secondary objectives were addressed: updating field inventory data collected in three different dates to the point clouds acquisition dates using single-tree growth models; comparing five selection methods and five regression methods in forest attributes modeling.

Study Area
The Aleppo pine (Pinus halepensis Mill.) forests under study are located in the Aragón region (Figure 1), Northeast Spain. This species represents 18.7% of the forested area, including semi-natural and reforested stands [48] and is well adapted to Mediterranean environmental conditions. Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 31

Study Area
The Aleppo pine (Pinus halepensis Mill.) forests under study are located in the Aragón region ( Figure 1), Northeast Spain. This species represents 18.7% of the forested area, including semi-natural and reforested stands [48] and is well adapted to Mediterranean environmental conditions. Study area with the location of forest inventory plots. High spatial resolution orthophotography from Spanish National Plan for Aerial Orthophotography spatial data infrastructure (SDI)service [54] is used as a backdrop.
In this area, the annual precipitation ranges from 350 mm to 1000 mm [49].The average annual temperature is 14° C, with cold winters and hot and dry summers. Aleppo pine forests are characterized by a hilly topography, with elevations ranging from 300 to 1150 m above sea level and slopes from 0° to 39°. The lithology of the study area varies from Miocene carbonate and marl sediments to limestone platforms and Mesozoic and Eocene limestone.
Some pine stands are natural, but most stands were planted approximately forty to sixty years ago. The evergreen understory includes species such as Quercus ilex subsp. rotundifolia, Quercus coccifera, Juniperus oxycedrus, Buxus sempervirens, Juniperus phoenicea, Rosmarinus officinalis, and Thymus vulgaris, among many others. Reforested areas usually present a low presence of hardwood species In this area, the annual precipitation ranges from 350 mm to 1000 mm [49]. The average annual temperature is 14 • C, with cold winters and hot and dry summers. Aleppo pine forests are characterized by a hilly topography, with elevations ranging from 300 to 1150 m above sea level and slopes from 0 • to 39 • . The lithology of the study area varies from Miocene carbonate and marl sediments to limestone platforms and Mesozoic and Eocene limestone. Some pine stands are natural, but most stands were planted approximately forty to sixty years ago. The evergreen understory includes species such as Quercus ilex subsp. rotundifolia, Quercus coccifera, Juniperus oxycedrus, Buxus sempervirens, Juniperus phoenicea, Rosmarinus officinalis, and Thymus vulgaris, among many others. Reforested areas usually present a low presence of hardwood species and poor diversity [50,51], while natural stands are structurally complex with a developed and diverse understory [52].

Forest Inventory Data
Forest inventory data was acquired through 147 plots from three field campaigns performed during June to July 2013, July to September 2014, and April 2016 (Figure 1), from now on cited as first, second, and third campaign, respectively. The sampling data fulfilled the statistical requirements [53], considered the size of the study area, and the variability of the pine forest in terms of terrain slope, canopy height, and canopy cover (CC) [54].
Field data from the first campaign were acquired in 53 circular plots with a 15 m radius. A Leica VIVA®GS15 CS10 GNSS real-time kinematic global positioning system (GPS) was used to collect the center of each plot. Tree dbh, for those trees with a dbh larger than 7.5 cm, was measured at 1.3 m height using a diameter tape. Green crown height and height of up to 4 randomly selected sample trees within each plot were measured using a Suunto®hypsometer. Thus, diametric class was considered when selecting the sample trees.
Field data from the second campaign were collected in 43 circular plots with a 15 m radius. The center of the plots was referenced using the same Global Navigation Satellite System (GNSS) receiver as in the first campaign. Tree dbh was measured with the same criteria as the first campaign using a Haglöf Sweden®Mantax Precision Blue caliper. The green crown height and the total height of all trees in the plot were measured using a Haglöf Sweden®Vertex instrument.
Field data from the third campaign were acquired in 51 circular plots. The center of the designated circular plots was measured using a Trimble®GNSS receiver. Field plots with a 5.6 m (3 plots), 8.5 m (23 plots), 11.3 m (17 plots), and 14.10 m (7 plots) radius were collected. Tree dbh was measured at 1.3 m using the same caliper as in the second campaign. The green crown height and the height of up to the 6 nearest trees to the plot center were measured using a Haglöf Sweden®Vertex instrument. The sample was completed to achieve 100 dominant stems ha −1 considering those with larger dbh.
The height for those trees not measured in the field plots was predicted using a height-diameter model developed from the sampled trees from all the field plots of the third campaign. Normality of the residuals, homoscedasticity, and independence or no auto-correlation in the residuals were verified for the linear regression fitted model.

Inventories Updating and Stand Variable Computation
Field data measurements were updated to year 2011 and 2016, which correspond to each ALS flight, to avoid any temporal lag between ALS-metrics and stand-level variables. The PHRAGON-2017 individual tree model was applied [55] through the forest simulator platform Simanfor [56]. This model enables tree-level distance-independent simulation of the development of Aleppo pine afforestations in Aragón. It consists of a set of equations for diameter over bark growth, diameter under bark growth, diameter under bark-diameter over bark ratio, generalized height-diameter relationship, volume over bark (taper equation) and crown ratio. In addition, it presents a survival model for the coming 10-year period and a classification tree for the regeneration of species of the genus Pinus, Quercus, and Juniperus, also in the coming decade. Explanatory variables included those related to tree size (diameter at breast height, total height), stand density (basal area, Hart-Becking index), dominant trees (dominant height, dominant diameter), competition (BALMOD) [57] and site quality (site index). Site index and dominant height evolution were estimated using the site index curves developed for natural Aleppo pine forests in the central Ebro valley [58].
Thus, when projecting future stand variables, diameter growth and survival equations were applied to every single tree in each plot, while the site index curve was used to forecast the future stand dominant height (and hence estimate the total height of each surviving tree). To reconstruct stand structure in the past, we need to deploy the diameter under bark growth equation, as it permits the use of the current stand features to predict the past growth of a tree (backdating procedure). Therefore, past tree diameters over bark are estimated through the diameter under bark -diameter over bark ratio-, while past dominant height can be calculated with the same site index curves, as they are dynamic, age-independent functions. Once every diameter and dominant height are known, the rest of the stand variables can be directly computed.
Seven forest inventory attributes were calculated from field data for each plot: stand density (N); basal area (G); squared mean diameter (Dg); dominant diameter (Do); dominant height (Ho); timber volume over bark of stem (V); and total tree biomass (W) [37] (Table 1). Thus, Ho is the mean height of the 100 trees per ha with largest dbh; Do is the mean dbh of the 100 trees per ha with largest dbh. V is estimated through the taper equation included in the PHRAGON-2017 individual-tree model [55]. W is computed as the sum of aboveground and belowground tree biomass using the Aleppo pine allometric equations developed by Ruiz-Peinado et al. [59].

ALS Data and Processing
The ALS data were acquired in 2011 and 2016 by the Spanish National Plan for Aerial Orthophotography (PNOA) [60]. The respective acquisition specifications are shown in Table 2. The point clouds, delivered in 2 × 2 km tiles in LAS binary file format, were captured with up to four returns measured per pulse. The x, y, and z coordinates were provided in European Terrestrial Reference System (ETRS) 1989 Universal Transverse Mercator (UTM) Zone 30 N. After the noise removal from the point clouds, a verification of the overlapping returns was performed considering vertical and horizontal displacements. Thus, overlapping returns were removed from 105 tiles for the 2011 ALS flight. The subsequent steps were evenly applied for both ALS campaigns. The multiscale curvature classification algorithm [61], implemented in MCC 2.1 command line tool, was used to classify ground points according to Montealegre et al. [62]. A digital elevation model (DEM) with a 1-m size grid was generated using the Point-Triangulated Irregular Network-Raster interpolation method [61], implemented in ArcGIS 10.5 software. This DEM was used for point height normalization. The point clouds were clipped to the spatial extent of each field plot. Then, a full suite of statistical metrics related to height distribution and canopy cover was calculated [63] using FUSION LDV 3.60 [64] software. A threshold value of 2 m height was applied to remove ground and understory laser hits before generating the ALS-derived variables according to Nilsson [65] and Naesset and Økland [66]. For better understanding of the results, the ALS metrics were classified into three macro classes and seven classes (see Table A1 in Appendix A). Canopy height metrics (CHMs) were subdivided into lower, mean, and higher height variables; canopy height variability metrics (CHVMs) were subdivided into variability and variability metrics derived from the L moments [47]; and canopy density metrics (CDMs) were subdivided into percentage of first or all returns, canopy relief ratio (CRR), and the ratio between all returns and total returns.  variability metrics (CHVMs) were subdivided into variability and variability metrics derived from the L moments [47]; and canopy density metrics (CDMs) were subdivided into percentage of first or all returns, canopy relief ratio (CRR), and the ratio between all returns and total returns.  Forest stand attributes modeling using the indirect approach was performed by two steps: (i) selection of the suitable ALS metrics using five selection approaches, and (ii) estimation of each stand attribute using five types of regression methods for 2011 and 2016 years ( Figure 2). Thus, each of the computed models have associated a selection approach, which determined the ALS metrics to be included in the models. As described by Domingo et al. [36], different selection methods were applied to choose the ALSmetrics that present the best relationship with the forest inventory attribute at plot-level: (i) Spearman's rank correlation coefficient considering a minimum positive and negative rho value of 0.5; (ii) stepwise selection using backward, forward, and bidirectional approaches; (iii) principal component analysis (PCA) using varimax rotation to better interpret the results [67,68]; (iv) last absolute shrinkage and selection operator (LASSO) [69]; and (v) all subset selection (ASS) implementing exhaustive, backward, forward, and sequential replacement (Seqrep) approaches [70].
The selection of ALS-metrics was restricted to a combination of up to four independent variables using the mentioned selection methods to avoid variable multicollinearity, overfitting [71], and within the purpose of generating parsimonious models for forest management.

Variable Selection and Attributes Modeling Using the Indirect Approach
Forest stand attributes modeling using the indirect approach was performed by two steps: (i) selection of the suitable ALS metrics using five selection approaches, and (ii) estimation of each stand attribute using five types of regression methods for 2011 and 2016 years ( Figure 2). Thus, each of the computed models have associated a selection approach, which determined the ALS metrics to be included in the models.
As described by Domingo et al. [36], different selection methods were applied to choose the ALS-metrics that present the best relationship with the forest inventory attribute at plot-level: (i) Spearman's rank correlation coefficient considering a minimum positive and negative rho value of 0.5; (ii) stepwise selection using backward, forward, and bidirectional approaches; (iii) principal component analysis (PCA) using varimax rotation to better interpret the results [67,68]; (iv) last absolute shrinkage and selection operator (LASSO) [69]; and (v) all subset selection (ASS) implementing exhaustive, backward, forward, and sequential replacement (Seqrep) approaches [70].
The selection of ALS-metrics was restricted to a combination of up to four independent variables using the mentioned selection methods to avoid variable multicollinearity, overfitting [71], and within the purpose of generating parsimonious models for forest management.
The estimation of forest stand attributes using an area-based approach and ALS data is usually done using either parametric (i.e. multiple linear regression) or non-parametric approaches such as regression trees, random forest, support vector machine, k-nearest neighbor, artificial neural network among others [72]. Five different regression methods, as described by Domingo et al. [36], were compared to estimate the seven forest inventory attributes (Table 1): multiple linear regression model (MLR), support vector machine (SVM), random forest (RF), locally weighted linear regression (LWLR), and linear model with a minimum length principle (MDL).
In the case of multiple linear regression model (MLR) normality, homoscedasticity, independence, and no auto-correlation of the residuals were verified. Logarithmic transformation of the dependent variables and the independent ones was also performed in those cases where statistical assumptions of linear regression were not fulfilled [73][74][75] or to improve the goodness-of-fit of the models. Support vector machine was computed using two kernel variants, linear (SVMl), and radial (SVMr) ones. Cost and gamma SVM parameters were tuned applying an interval of 1-1000 and 0.01-1, respectively. Random forest was implemented in R using "randomForest" [76] and "caret" packages [77], including "corr.bias" parameter to minimize bias effects. The model was tuned by applying a number of trees to growth (ntrees) within the interval 1-3000 and a number of variables to divide the nodes between 1 and 3. Locally weighted linear regression and MDL tree structures were computed using the R package "CORElearn" considering up to four ALS metrics. Model computation required the splitting of the original sample into a training set of 75% of the cases (110 plots) and a test set of 25% of the cases (37 plots). Validation was performed for all the models, being executed 100 times for those methods with associated randomness as SVM, RF, LWLR, and MDL to increase robustness in the results [78]. Furthermore, data were normalized in values ranging from 0 to 1 in order to avoid weights saturation according to Görgens [79].
Statistic performance of each computed model was reported including root mean square error (RMSE), relative RMSE respect to the mean (%RMSE) and bias. Differences between models were determined by using the Friedman nonparametric test according to the RMSE of each plot [80]. Furthermore, the Nemenyi post-hoc test was applied to determine whether differences were statistically significant, with a significance level of 0.05 [81]. This test was applied only when the null-hypothesis of the Friedman test was rejected, thus implying non-equivalence between models.

Assessment of Temporal Transferability by Applying a Direct Approach
The temporal transferability of models were assessed by three steps: (i) selection of the best regression model previously generated by the indirect approach for each forest stand attribute and year (2011 and 2016); (ii) extrapolation of the selected models from 2011 to 2016 ALS data, using the same variables and model parameters, and inversely from 2016 to 2011 by using the direct approach; (iii) performance comparison between extrapolated models for both years ( Figure 2). Thus, Friedman and Nemenyi tests were applied for selecting the best regression model for each year (step i) and for selecting the best transferable models for both years (step iii).

Field Plot Computation
Equation (1) was used for estimating tree height (ht) for those trees not measured in the field plots for the first and third campaign. Model performance for the ht model was as follow: RMSE of 0.80 m and R 2 of 0.93.
where ht is tree height (m), dbh is the diameter at breast height (cm), H 0 and D 0 are the dominant height and dominant diameter, as defined in Section 2.3. Table 3 shows a summary of the forest inventory attributes obtained from the field plot data. The N values of the inventoried plots ranges from 99.03 to 3200 stems ha −1 and G ranges from 0.82 to 58.89 m 2 ha −1 , presenting also a variety of diameters, from 9.21 to 47.96 cm. V and W data show also a high range of values with a standard deviation in both cases higher than 60 tons ha −1 . The high range and standard deviation values of the forest inventory attributes show the variability that characterizes Aleppo pine forest in Aragón region.  Table 4 shows a summary of the estimated field plot attributes using single-tree growth models for each measured stand variable and both years. N shows a general decrease in the number of stems ha −1 , which may be caused by tree growth, resulting in an average of N change of 10.

Variable Selection
This section includes the results of the selection of ALS variables for the seven estimated forestry attributes modeled with the different regression methods. Figure 3A synthetizes the performance of the analyzed selection methods for each forest stand attributes by year. All subsets regression Seqrep (ASSs) was the most powerful selection method. Spearman's rank correlation (rho) coefficient also showed good results, especially for selecting N, G, Dg, and W in 2011. All subsets regression Exhaustive (ASSe) and Stepwise (Step b&f) were both good selection methods for estimating G, Ho, and V. However, lasso selection (LASSO), all subsets regression Backward (ASSb) and all subsets regression Forward (ASSf) have been less utilized. The stepwise forward and PCA selection methods have not been included in Figure 3 as they did not determine the best variables for modeling in any of the cases. For detailed information of the selection methods performance, see Tables A2-A15 in the Appendix A.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 31 Figure 4 shows the ALS selected metrics for estimating the forest inventory attributes for both 2011 and 2016 years. As mentioned in Section 2.4, for better understanding of results, the ALS metrics are classified into groups (see Table A1 in Appendix). In general, higher height variables, variability, and variables related to the ratio between all returns and total returns were included in most of the models, while height variables and variability L moments were less demanded. Comparing the different estimated forest attributes, some differences can be observed. N, Dg, and Do estimations included higher height variables, variability metrics, and CDM metrics. Ho estimation usually required higher height variables and variability metrics, while CDM metrics were not included. G, V, and W estimations included either lower or higher height metrics, variability metrics, and CDM ones.  Figure 3A was six for all the stand attributes except for volume (V) and total tree biomass (W), which have a maximum number of five models. Maximum number of models in Figure 3B is 14, seven for each year and stand attribute, for all regression methods. Rho stands for Spearman Rank; ASSs stands for All Subsect Selection Seqrep; ASSe stands for All Subsect Selection Exhaustive; ASSf stands for All Subsect Selection Forward; ASSb stands for All Subsect Selection Backward; and Step. b&f stands for Stepwise Selection Both Backward and Forward.  Figure 3A was six for all the stand attributes except for volume (V) and total tree biomass (W), which have a maximum number of five models. Maximum number of models in Figure 3B is 14, seven for each year and stand attribute, for all regression methods. Rho stands for Spearman Rank; ASSs stands for All Subsect Selection Seqrep; ASSe stands for All Subsect Selection Exhaustive; ASSf stands for All Subsect Selection Forward; ASSb stands for All Subsect Selection Backward; and Step. b&f stands for Stepwise Selection Both Backward and Forward. Figure 3B depicts the performance of selection methods associated to each regression method without considering the year of the model. The ASSs was the most powerful method to select the best ALS metrics when using the MDL, LWLR, SVMr, and SVMl regression methods. Furthermore, rho coefficient was the most powerful one when using the MLR and RF regression methods. The ASSe and Step b&f both were also good selection techniques for almost all the regression methods, excluding MLR. However, LASSO, ASSb, and ASSf were less effective. Figure 4 shows the ALS selected metrics for estimating the forest inventory attributes for both 2011 and 2016 years. As mentioned in Section 2.4, for better understanding of results, the ALS metrics are classified into groups (see Table A1 in Appendix A). In general, higher height variables, variability, and variables related to the ratio between all returns and total returns were included in most of the models, while height variables and variability L moments were less demanded. Comparing the different estimated forest attributes, some differences can be observed. N, Dg, and Do estimations included higher height variables, variability metrics, and CDM metrics. Ho estimation usually required higher height variables and variability metrics, while CDM metrics were not included. G, V, and W estimations included either lower or higher height metrics, variability metrics, and CDM ones.

Indirect Approach
The regression models to estimate Dg, Do, G, Ho, N, V, and W for 2011 and 2016 years using the indirect approach are summarized in Tables A2-A15 of the Appendix. Figure 5 summarizes the %RMSE respect to the mean of the different regression methods for estimating the forest inventory attributes. Models developed for 2016 ( Figure 5B) present generally higher accuracy than the ones generated for 2011 ( Figure 5A). The point density of ALS datasets may determine these differences in accuracy. The SVMr shows the lower RMSE when modeling all the analyzed stand attributes in 2011 ( Figure 5A). In this year, SVMl was the second best model when estimating Do, G, Ho, and W; MDL when estimating N and V; and RF when estimating Dg. The MLR regression method was not computed for V and W, as statistical assumptions of linear regression were not fulfilled, even considering logarithmic transformation. The MLR shows the lowest accuracy when estimating Do, G, Ho, and N.

Indirect Approach
The regression models to estimate Dg, Do, G, Ho, N, V, and W for 2011 and 2016 years using the indirect approach are summarized in Tables A2-A15 of the Appendix A. Figure 5 summarizes the %RMSE respect to the mean of the different regression methods for estimating the forest inventory attributes. Models developed for 2016 ( Figure 5B) present generally higher accuracy than the ones generated for 2011 ( Figure 5A). The point density of ALS datasets may determine these differences in accuracy. The SVMr shows the lower RMSE when modeling all the analyzed stand attributes in 2011 ( Figure 5A). In this year, SVMl was the second best model when estimating Do, G, Ho, and W; MDL when estimating N and V; and RF when estimating Dg. The MLR regression method was not computed for V and W, as statistical assumptions of linear regression were not fulfilled, even considering logarithmic transformation. The MLR shows the lowest accuracy when estimating Do, G, Ho, and N. The SVMr is the best model for estimating N, V, W, Dg, and G in 2016 ( Figure 5B). However, in this year, SVMl and MDL outperformed SVMr when estimating Do and Ho, respectively. In 2016, MLR shows the lowest accuracy when estimating G, Ho, and N.
Friedman tests shows that the models are not equivalent, with a p-value lower than 0.05 when testing whether there were statistically significant differences between regression methods for 2011 and 2016 years. Thereby, the post-hoc Nemenyi test indicates that no statistically significant differences exist between the methods, with 95% of probability, except for MLR. In this sense, statistically significant differences were found when comparing models in the following cases: between MLR and SVMr models when estimating Do and G for 2011; between MLR and MDL models when estimating Ho for 2011; between MLR and MDL models when estimating Dg for 2016; and between MLR and all the generated models when estimating G, Ho, and N for 2016.

4. Direct Approach
The SVMr was established as the regression method for analyzing how models fitted at 2011 perform at 2016, and inversely, following a direct approach to analyze temporal transferability. This method resulted the best estimator for all the models generated in 2011 and for the majority of forest attributes modeled in 2016. Table 5 summarizes the best-selected SVMr models fitted in 2011 and the ones extrapolated to 2016 by using the same ALS metrics and model parameters. Table 6 summarizes the best-selected SVMr models fitted in 2016 and extrapolated to 2011 by using the same ALS metrics and model parameters. Furthermore, scatter plots of observed and predicted values for the analyzed forest stand attributes for both years are included in Figures A1 and A2 of Appendix.
In the case of the models fitted in 2011 and extrapolated to 2016 (Table 5), the %RMSE after validation ranges from 8.54 to 42.43% and R 2 ranges from 0.64 to 0.93 within the different stand attributes. As it is shown in Table 5, models are transferable. In fact, the average %RMSE differences between the fitted and the extrapolated models is 3.87%. Dg, Do, Ho, V, and W estimations for 2011 models have higher %RMSE than the one for models extrapolated to 2016. However, N and G models show higher %RMSE for the 2016 extrapolated ones.
In the case of the models fitted in 2016 and extrapolated to 2011 (Table 6), the %RMSE in the validation sample ranges from 8.83 to 48.09% and R 2 ranges from 0.55 to 0.92 within the different stand attributes. These models also show good temporal transferability, being the average %RMSE differences between the fitted and the extrapolated model 5.85%, even lower than the models fitted The SVMr is the best model for estimating N, V, W, Dg, and G in 2016 ( Figure 5B). However, in this year, SVMl and MDL outperformed SVMr when estimating Do and Ho, respectively. In 2016, MLR shows the lowest accuracy when estimating G, Ho, and N.
Friedman tests shows that the models are not equivalent, with a p-value lower than 0.05 when testing whether there were statistically significant differences between regression methods for 2011 and 2016 years. Thereby, the post-hoc Nemenyi test indicates that no statistically significant differences exist between the methods, with 95% of probability, except for MLR. In this sense, statistically significant differences were found when comparing models in the following cases: between MLR and SVMr models when estimating Do and G for 2011; between MLR and MDL models when estimating Ho for 2011; between MLR and MDL models when estimating Dg for 2016; and between MLR and all the generated models when estimating G, Ho, and N for 2016.

Direct Approach
The SVMr was established as the regression method for analyzing how models fitted at 2011 perform at 2016, and inversely, following a direct approach to analyze temporal transferability. This method resulted the best estimator for all the models generated in 2011 and for the majority of forest attributes modeled in 2016. Table 5 summarizes the best-selected SVMr models fitted in 2011 and the ones extrapolated to 2016 by using the same ALS metrics and model parameters. Table 6 summarizes the best-selected SVMr models fitted in 2016 and extrapolated to 2011 by using the same ALS metrics and model parameters. Furthermore, scatter plots of observed and predicted values for the analyzed forest stand attributes for both years are included in Figures A1 and A2 of Appendix A. In the case of the models fitted in 2011 and extrapolated to 2016 (Table 5), the %RMSE after validation ranges from 8.54 to 42.43% and R 2 ranges from 0.64 to 0.93 within the different stand attributes. As it is shown in Table 5, models are transferable. In fact, the average %RMSE differences between the fitted and the extrapolated models is 3.87%. Dg, Do, Ho, V, and W estimations for 2011 models have higher %RMSE than the one for models extrapolated to 2016. However, N and G models show higher %RMSE for the 2016 extrapolated ones.
In the case of the models fitted in 2016 and extrapolated to 2011 (Table 6), the %RMSE in the validation sample ranges from 8.83 to 48.09% and R 2 ranges from 0.55 to 0.92 within the different stand attributes. These models also show good temporal transferability, being the average %RMSE differences between the fitted and the extrapolated model 5.85%, even lower than the models fitted in 2011 and extrapolated to 2016. All the models fitted in 2016 for the analyzed stand attributes present lower RMSE than the ones extrapolated to 2011.
The comparison between fitted models generated for either 2011 or 2016 and the extrapolated ones were assessed using Friedman and post-hoc Nemenyi tests. Friedman test shows that the models for the analyzed stand attributes are not equivalent with a p-value lower than 0.05 when testing whether there were differences: (i) between models fitted in 2011 and the ones extrapolated to 2016; (ii) between models fitted in 2016 and the ones extrapolated to 2011; (iii) between models fitted in 2011 and models fitted in 2016; (iv) between models fitted in 2011 and models extrapolated to 2011; (v) between models fitted in 2016 and models extrapolated to 2016. Thereby, the required application of the post-hoc Nemenyi test indicates that no statistically significant differences exist between the methods for the proposed hypothesis, with 95% of probability.

Discussion
Airborne laser scanning is considered the best technology for mapping 3D vegetation structures [3] allowing the measurement of fine-scale forestry metrics [82]. Multi-temporal ALS data has been less explored as the availability of two or more LiDAR surveys in the same area is still limited. Nevertheless, several studies have used multi-temporal small-footprint ALS to estimate total tree biomass or carbon dynamics [3,9,[21][22][23][24][25][26][27], volume changes [17,21], height growth [14][15][16][17][18][19], and basal area [17]. This study estimates seven forest attributes (N, G, Dg, Do, Ho, V, W) using bi-temporal low-point density ALS data in Mediterranean Aleppo pine heterogeneous forests. The high number of field plots has allowed estimating the seven mentioned forest attributes for 18.7% of the forested area in Aragón, providing results at a regional scale. Moreover, model temporal transferability was demonstrated which could improve forest management in a cost-effective way in Mediterranean Aleppo pine forests.
Multi-temporal LiDAR estimations of forest attributes requires the support of accompanying field surveys [3] being desirable to have them corresponding to LiDAR surveys [9]. Field surveys are cost and time demanding specially when acquiring a high number of plots. The use of specific individual-tree growth equations, derived from dbh growth by extracting tree cores or from interval or permanent plots, is a good way to get value from field plot inventories acquired between different ALS surveys. Diameter at breast height and height growth values from general yield tables from the Spanish National Forest Inventory have been applied in other studies for estimating total tree biomass in Aleppo pine forests [36]. Nevertheless, in this work specific single-tree growth models, generalized height-diameter curves and taper equations were used to estimate all stand attributes in the measured field plots at two different points in time, which produces more accurate results, particularly when predicting at short term [83].
The use of selection methods reduces variable collinearity, modeling time and increases model parsimony. All subset selection Seqrep was the most powerful selection method in agreement with Hansen et al. [84] who also used similar best subset regression procedures to estimate biomass with ALS data. Spearman rho coefficients, proposed as a good tool for determining the relationships between ALS and field metrics by Kristensen et al. [85], also showed a good result, agreeing with our previous studies [37]. Furthermore, our results agree with García-Gutiérrez et al. [78], who found that stepwise was a powerless technique. Accordingly, the use of automatic selection methods such as ASSs is proposed when using MDL, LWLR, SVMr, and SVMl regression methods in Mediterranean Aleppo pine forest. Nevertheless, comparison between selection methods should be considered when working with other forest types or species. In this sense, Rho coefficients should be considered specially when using MLR and RF regression methods and PCA should be taken into account for a first attempt to reduce collinearity as proposed by Naesset [54] and Cao et al. [9], but afterwards another selection approach should be considered before modeling.
The most selected types of ALS metrics for estimating the seven analyzed stand attributes were higher height variables, variability ones and the ratio between all returns and total returns, while dominant height and variability L moment variables where less demanded. Ho estimation usually required the inclusion of high height percentiles as concluded by Naesset and Gobakken [17]. V and W estimations normally included either lower or higher height variables, variability metrics, and/or CDM ones as proposed by Silva et al. [86] and Hopkinson et al. [87].
The comparison between regression methods shows that SVMr had the lowest RMSE when estimating the majority of the analyzed stand attributes for both dates, except for Do and Ho when using 2016 data. These results match with García-Gutiérrez et al. [78] and Guerra-Hernández et al. [38,88], which obtained higher accuracies when using non-parametric regression methods. Different results were found in our previous studies [36,89] as MLR slightly outperformed other nonparametric methods when estimating total tree biomass in Aleppo pine forests, but no statistically significant differences were found. Thus, in agreement with Gagliasso et al. [90], a high number of field plots may have boosted machine-learning performance. Furthermore, the broad range and standard deviation values of the field plot data that characterizes Aleppo pine forest at a regional scale is notoriously higher than in our previous studies. Thus, although logarithmic transformation of the dependent and independent variables was carried out, most of the data was not normally distributed. The limitation on using the best-suited ALS metrics, as most of them were not normally distributed, generates a considerable decrease of accuracy in MLR model performance. In this sense, comparison between regression methods is desirable, especially when working with big datasets in heterogeneous forest stands as the Mediterranean ones.
The comparison of direct and indirect approaches allowed us to assess model temporal transferability between 2011 and 2016. The direct approach was computed when extrapolating 2011 models to 2016 and inversely. The indirect approach has shown slightly better results when estimating N, G, Dg, Do, and V. Direct approach showed slightly better results in the estimation of W when extrapolating 2016 model to 2011, but not inversely. Similar results were found for the estimation of Ho when extrapolating 2011 model to 2016 data, but not inversely. Comparisons with previous studies cannot be done for N, Dg, and Do, as these attributes have not been previously estimated using multi-temporal data. Comparisons between Ho, G, and V results are neither possible as Naesset and Gobakken [17] performed only the indirect approach. Regarding W estimations, several results have been obtained using different regression methods, and even in our work both direct and indirect approaches performed in a different way when extrapolating the first-year model (2011) or the second one (2016). The indirect approach obtained better results when estimating W for 2011 data agreeing with Zhao et al. [3] and Meyer et al. [25]. Our results also agree with Cao et al. [9], Skowronski et al. [24], and Bollandsås et al. [26], which detected slightly better performance of direct approach. Similar results have been found in our study when extrapolating 2016 model to 2011. In general, the SVMr regression method shows good temporal transferability between both ALS acquisitions with average %RMSE differences for all the modeled stand attributes of 2.13% for 2011 and 1.58% for 2016.
Models generated using 2016 data (1.25 points m −2 ) showed generally higher accuracy than 2011 ones (0.64 points m −2 ). However, no statistically significant differences were found between the best-fitted models for each year. In agreement with Cao et al. [9], point density may influence model performance but did not significantly affect the estimations of forest attributes as point clouds has a consistent vertical pattern. According to Hudak et al. [27], the relatively large size of the sample plots is considered sufficient for generating canopy height metrics. Thus, the results confirm, as other previous studies based on low-density ALS from the Spanish National Plan for Aerial Orthophotography data (i.e., [33,[36][37][38]), that this information is an accurate and economic alternative to perform forest inventories when higher point density data are not available.
Overall, the use of low-point ALS data for two dates and single-tree growth models for generating temporally-concomitant field data provides accurate estimations of forest stand attributes in Mediterranean Aleppo pine forests. The indirect approach produced higher precision, but the direct approach, within those conditions, may reduce fieldwork and time of model parametrization. When using a direct approach it would not be necessary to create one model for two different points in time, as it will be possible to extrapolate a model generated for one date (validated with field data) to another date. Furthermore, the number of revisited field plots can be dismissed or even not required, when the time between ALS surveys is not large [28]. This will benefit not only forest managers but also enterprises devoted to forest inventories. The use of direct method and the possibility of model temporal transferability generates new alternatives to calibrate future ALS captures with a lower number of field plots and helps in designing the temporal gap between flights. Single-tree growth models constitute a useful and robust alternative to update field data to a point in time, allowing to accurately estimate forest inventory parameters with the use of ALS data. Future research using multi-temporal ALS data should focus on the inclusion of inference models to better understand uncertainties as well as on the analysis of field plot size and saturation effects in model accuracy. Furthermore, the analysis of forests structural biodiversity changes caused by wildfires or the fusion of ALS data with multi-temporal passive remote sensing series or unmanned aerial vehicle (UAV) point clouds may help to monitor forest dynamics over time.

Conclusions
Multi-temporal ALS data may improve forest management and planning, providing accurate forest inventory attribute estimations for different points in time. The results illustrate the usefulness of bitemporal low-point density ALS data and single-tree growth models, when lacking temporallyconcomitant field campaigns, to accurately estimate seven forestry attributes, using an area-based approach. All subsets regression Seqrep was the most powerful selection method, followed by rho coefficient, to generate parsimonious models. Higher height metrics, canopy height variability, and canopy density variables were the most selected ALS-metrics, while mean height variables and variability L moments were less demanded. The SVM with radial kernel outperformed the analyzed non-parametric and multivariate linear regression methods for estimating all forest inventory attributes except from Do and Ho when using 2016 data. Thus, machine-learning performance may have been boosted by forest heterogeneity and an elevated number of field plots.
This study has assessed model temporal transferability by comparing direct and indirect approaches for the estimation of seven forestry attributes. Indirect approach have produced slightly more accurate results than the direct approach, but average %RMSE differences between both approaches for all modeled stand attributes ranged from 2.13% in 2011 to 1.58% in 2016. Thus, mixing the direct approach with single-tree growth methods provides a suitable alternative to reduce fieldwork and enhance ALS technology as a good tool for estimating forest attributes in two different dates. The utility of multi-temporal ALS data and the combination with multi-temporal series from passive remote sensing and UAV point clouds derived by using photogrammetric techniques would have great value for forest management and planning.  We thank the Spanish National Geographic Information Centre (CNIG) and the Geographic Institute of Aragón (IGEAR) for providing the ALS data. We are also grateful to the material resources provided by the Centro Universitario de la Defensa de Zaragoza (CUD). We express our appreciation to Alberto García-Martín, Jesús Cabrera, and Francisco Escribano for their help during the fieldwork.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Summary of the airborne laser scanning (ALS) computed metrics including the abbreviations, classes, and macro-classes defined.

Macro-Classes Classes ALS Computed Metrics Abbreviations
Canopy height metrics (CHM) percentage of first returns above the 2.00 % first ret. above 2.00 percentage of all returns above the 2.00 % all ret. above 2.00 percentage of first returns above the mean % first ret. above mean percentage of first returns above the mode % first ret. above mode percentage of all returns above the mean % all ret. above mean percentage of all returns above the mode % all ret. above mode Canopy relief ratio CRR All returns Total returns-1 All returns above 2.00 divided by the total first returns × 100 (All ret. above 2.00)/(total first ret.) × 100 All returns above mean divided by the total first returns × 100 (All ret. above mean)/(total first ret.) × 100 All returns above mode divided by the total first returns × 100 (All ret. above mode)/(total first ret.) × 100     Table A13. Summary of the V models using 2016 ALS data. Validation results in terms of RMSE (m 3 ha −1 ), %RMSE, and bias (m 3 ha −1 ) and R 2 . SM refers to selection method; Step. stands for Stepwise both and forward; SVMr. refers to support vector machine with radial kernel; SVM l. refers to support vector machine with linear kernel; ret. refers to returns.   Figure A1. Scatterplot of predicted vs. observed values of the forest stand variables using the bestselected SVM with radial kernel 2011 models and 2016 extrapolated models. Dg and Do are expressed in cm; G is expressed in m 2 ha -1 ; H is expressed in m; N is expressed in stems ha -1 ; V is expressed in m³ ha -1 ; W is expressed in tons ha -1 .