Assessing understory development in forest plantations using laser imaging detection and ranging ( LiDAR ) Evaluación del desarrollo del sotobosque en plantaciones forestales mediante LiDAR

Forestry plantations are expected to be managed in ways to conserve biodiversity while producing goods and services. This goal implies a signifi cant challenge as plantations tend to reduce species richness. The presence of well developed understory enhances the value of plantations as habitat for native fauna. Here, we develop a straightforward method to assess the availability of understory in forestry stands using laser imaging detection and ranging (LiDAR) data and aerial RGB high resolution images. Based on fi eld and airborne acquired data for Pinus radiata stands in central Chile, the digital crown model (DCM), derived from the subtraction of the digital terrain model (DTM) from the digital surface model (DSM) is a more reliable predictor of understory height that variables like terrain slope, aspect, plantation age and canopy height in forests and plantations which have not complete closed canopy. The correlation between DCM and understory though decreases while the actual height of the plantation canopy increases, rendering DCM a conservative estimate of understory development. The use of DCM will allow a fast and cost/effective estimate of habitat suitability in forestry plantations.


INTRODUCTION
Forestr y plantations covers over 265 million ha worldwide, growing 2 % yearly (FAO 2011).Landscape transformation brought about by these plantations has impinged upon biological diversity, particularly planted forests based on introduced species (Carnus et al. 2006).Indeed, the Convention of Biological Diversity (CDB) has set as a decennial goal to enhance the role of forestr y plantations in the conser vation of biodiversity (CBD 2010).Increasing evidence suggests that structurally complex plantations, with a well developed understory could be a way forward to achieve this goal.In fact, species richness and abundance tend to be higher in planted forests with well developed understor y than those devoid of such vegetation, therefore minimizing the reduction of species richness of soil insects, herpetozoans, mammals and birds among other taxa (see Simonetti et al. 2012 for a review).
In this framework, to accurately and cost/ effectively assess the presence of understor y vegetation in planted forests is mandator y.Precise assessment of understor y coverage including its spatial heterogeneity enables quantification and characterization of the amount of suitable habitat that planted forests could provide, contributing to mitigate the impacts of landscape transformation.
Chilean conifer plantations are a case in point.Cur rently, 1.5 million ha are a monoculture of Pinus radiata D. Don.Stands with well developed understor y hold more species of insects, mammals and birds than those with scarce or no undergrowth.Among species using plantations, at least temporarily, are endangered species like Leopardus guigna (Molina), a fact that contributes to achieve CBD' goals (Simonetti 2006, Estades et al. 2012).Therefore, assessing how much of the area covered by plantations could be suitable as surrogate habitat for native fauna would provide solid managerial information to properly engage the forestry industry in biological conservation.
Here we develop a method to assess understory development using laser imaging detection and ran ging (LiDAR) data and aerial RGB high resolution images (see Baltsavias 1999, Wher & Lohr 1999and Elmqvist et al. 2001 for details regarding the use of these techniques).
LiDAR data, in conjunction with various sources of forest inventories data, can be successfully used to quantify different aspects of forests 3-D structure, such as volume and biomass estimates and canopy height profi les, among other features (Koch & Dees 2008).Hence, small-footprint LiDAR may be used to quantify components of forest structure needed to assess animal-habitat relationships (see review by Vierling et al. 2008).In fact, airborne laser scanning data can be used to predict habitat quality and to map species distributions as a function of habitat structure (Bradbur y et al. 2005, Goetz et al. 2007).Regarding the assessment of understor y vegetation, LiDAR data can be utilized for to derive both a variety of environmental factors which account for the presence of understor y vegetation, including canopy structure, forest successional stage and topography (e.g., Falkowski et al. 2009) and, to characterize height and cover of the understory vegetation either as woody vegetation or suppressed trees integrating only leaf-on and leaf-off LiDAR (Hill & Broughton 2009).
One critical step in using point clouds to derive digital terrain and surface models (DTM and DSM, respectively) is the identification and classification of ground returns, which impinges upon the quality and precision of any LiDAR product.Errors in the estimation of the ground level will affect the subsequent estimation of an object above it, including all types of vegetation.Many algorithms have been developed to generate DTMs using LiDAR data.These algorithms can be classifi ed into two categories based on type of data used: (a) point clouds and (b) raster range image.A DTM can be constructed from the ground points, and a digital surface model (DSM) can be derived from the highest points within a defi ned grid box (Hollaus et al. 2006).An object height model can be calculated then by subtracting DTM from the DSM.In a forest environment, without human built structures, all objects are assumed vegetation given a digital crown model (DCM).In this study, we deal with the use of DTM and DSM to estimate understory height in forestr y plantations.Our aim is to develop a simple step algorithm for processing these types of data in order to advance a simple novel way to estimate understory development, moving from comparisons of observed versus estimated values at plot level to the assessment of obser ved versus estimated understor y heights surfaces, which we obtained applying geo-statistical interpolation on fi eld and LiDAR data, respectively.This procedure will enable the assessment of the availability of suitable habitat for native fauna in forestry plantations, which could be used temporarily used for the native fauna, contributing to the fulfi llment of the CBD targets.

Study area
The

Field data collection and processing
Field data was gathered through systematic sampling of a grid of clusters spaced 200 m x 200 m over all types of vegetation.Each cluster contains fi ve subplots, one central and four neighboring ones 30 m apart, placed in the four geographic directions (N, S, W and E; Fig. 1).The coordinates of the plot center were located using a Garmin double frequency GPS Topcon GRS-1 (WGS84 UTM18s).We measured DBH, height, and crown ratio of all trees above 5, 10 and 20 cm of DBH in concentric rings of 2, 4 and 8 m in every plot.Understory was assessed only at the 4 m radii ring, recording species, average height, mean coverage (%) and plant density.A total of 98 clusters, comprising 490 plots were sampled.Using this data, a raster was obtained by applying ordinary kriging, representing understory height (H), as an estimate of its development.

LiDAR and image data acquisition
We decided to use LiDAR because the management of pine plantations, including pruning and thinning, determines that the canopy is not completely closed, which in turn allows both enough light to reach the ground understory development and returns of LiDAR pulses.
Discrete waveform laser scan data were acquired in March 31, 2011 by Digimapas Chile Ltd. (http://www.digimapas.cl)using a set Harrier 56/G4 Dual System mounted on a Piper PA-24 Comanche.A digital image acquisition sensor (VIS) and LiDAR scanning were set together and the fl ight was conducted with a nominal height over the ground of 580 m AGL without GPS errors, with an average speed between 180 to 210 km h -1 and leading to an average point density of 4.64 points per square meter (see Table 1 for fl ight and system details).Additionally, we acquired a aerial image (VIS), with the three channels (red, green and blue) in the visible range of the electromagnetic spectrum, with 1 m spatial resolution, obtained at the same time of LiDAR data and therefore it was truly orthorectifi ed.

LiDAR and image data processing
To process LiDAR point clouds into XYZ Ascii format, we used LasTools (Isenburg & Shewchuk 2010).First, we convert LiDAR data into *.las format and then we used the multiscale curvature algorithm (mcc-lidar) to classify ground points only (Evans & Hudak 2007).We used 1.5 and 0.3 as scale and threshold parameters.The optimal scale parameter depends on, fi rst, the scale of the objects on the ground (i.e.rocks, trees) and second, on the sampling interval (post spacing) of the LiDAR data.LiDAR sensors are capable of collecting high density data (e.g., 8-10 pulses m -2 ) which translate to 0.35 m pulse.This is 1/sqrt (8 pulses m -2 = 0.35 m pulse).Sparser LiDAR data (e.g., 0.25 pulses m -2 ) translate to a spacing of 2.0 m pulse (1/sqrt(0.25pulses m -2 ) = 2.0 m pulse).Evans & Hudak (2007) recommended 0.3 as a good starting value to try for the curvature threshold in forest landscapes.Once the ground points were classifi ed, we used the interpolation to raster LasTools (las2dem.exe)to obtain a DTM, using only ground points, and a DSM using all points.Further processing, subtracting DTM from DSM allow us to obtain a DCM.Complementary, the VIS image was used to identify and classify all stands in the study area into fi ve stand types: pine plantations (53 %), eucalyptus plantations (6 %), native forest (29 %), mixed stands (6 %), and other types (6 %).A polygon shape fi le was fi nally obtained for all former types.The DCM was sliced into eight height stratum as follows (in meters): 0-0.5, 0.5-1.0,1.0-2.0,2.0-4.0,4.0-8.0,8.0-16.0 and 16.0-32.0.Each one was then transformed into a raster containing only pixels within the correspondent height class.These pixels were separated into two sets, 60 % for training and performing ordinary kriging and 40 % for validation.The interpolated rasters were named H050, H100, H200, H400, H800, H1600 and H3200 to refl ect the height stratum they represented.Slope (SLOPE) and aspect (ASPECT) were directly calculated from the DTM data.Finally, the year of plantation (YEAR) was added as an additional variable to complete a fi nal set of 11 predictors.These variables are the most frequently used to model forest and plantation attributes (see Hernández et al. 2007).

Modeling
To be able to model the relations between all preprocessed predicted variables and the target variable, understory height (H), we sampled all rasters using 2000 plots of 4 m radii, randomly located over the study area.This way we obtained a master table that contains 2000 records and the 11 predictors plus understory height (H).Using this dataset we analyzed correlation among variables and fi tted multiply linear models.

Interpolation of fi eld and LiDAR data
The interpolation of heights using LiDAR data available for the seven strata height classes generates not only prediction sur faces but also errors, i.e. uncertainty surfaces, giving an indication of how good the predictions are.The error for all height classes has low values, near to zero, and also constrained standard deviation except for higher stratum where their values increase.The same behavior was obser ved for the kriging standard error, suggesting estimates are accurate (Table 2).

Correlations between understory height and predictors
The pr edictor DCM had the str ongest correlation to understor y height, exhibiting values that decrease as stratum height increase (Table 3).There is also a signifi cant (α = 0.05) positive Pearson's product moment correlation coeffi cient (> 0.3) between the overall stands (All) and predictors DCM, SLOPE, H400, H1600 and H3200.The same tendency was observed for stand of 0-0.5 m dominant height and DCM, and for stand of 8.0-16.0m and 4.0-8.0m dominant height.

Linear model
Figure 2, bottom left and right, shows the relation between the residuals of estimates and slope.There is evidence of an additive effect on residues when high slopes and high heights of plantations are combined that can be seen in the middle zone of the study area.On the contrary, in the bigger native forest stand, in the south-east area, there are also high slopes but the behavior of residues are the opposite, concentrating low values, probably because of the open canopy structure of this vegetation.Aspect was codifi ed into eight classes following a clockwise fashion from 1 (N) to 8 (NW) and therefore, the overall negative correlation is af fected by this codification.An interested behavior can be seen for the year of plantation variable (YEAR), which tends to increase its negative correlation with H as the age of the plantation increases.When the dominant canopy develops into denser structure, the understory exhibit less development.
Table 4 illustrates the linear model summary for the general form H = f (predictors).The sampling points were selected for each model according to their dominant height stand, from DCM, as indicated, fi tting models for each one.Most predictors showed strong lineal relation according to the stand dominant height under consideration.When the actual dominant height is low, i.e. less than 2 m, height stratum predictors have low explanation power.The overall model, including all ten predictors gave the best results.Strongly significant linear relation was also found at 99 %, for model 0 - 0.5 in DCM predictor, and for model 8.0 -16.0 in SLOPE predictor.The associated residues showed a slightly positive bias that make estimates to have a tendency to produce lower values than observed ones (Fig. 3).Stratum height variables, from H050 to H3200, did not provide strong correlation in any of the stands.This fact could be related to an inadequate selection of height classes, which separate height values in arbitrar y classes.Nevertheless, when using them as predictors in linear models, most of them showed signifi cant statistical relation to understory height (Table 4).All models exhibit strong lineal relation for a subset of predictors but, different ones, according to the dominant height of the stand under consideration.When the actual dominant heights are low, i.e. less than 2 m, the height stratum predictors have very little explanation power, which is probably related to the absence of either understor y or plantation species development.Both components tend to be Cuadriculados interpolados de la altura media del sotoboques observada (arriba-izquierda) y predicha (arriba-derecha).Abajoizquierda: mapa de los residuos estandarizados de modelo lineal general.Abajo-derecha: mapa de clases de pendientes.Summary of linear models for the general form H = f (predictors).Every column includes a different set of predictors as height strata varies.Signifi cant correlation at α = 0.001 are denoted by '***', α = 0.01 by '**', α = 0.05 by '*'and α = 0.1 by '.'; -indicates it was nor included in the indicated model.
mixed up and their discrimination is diffi cult to achieve.The overall model, including all predictors gave the best results, and can be use in understory mapping.

DISCUSSION
The Convention on Biological Diversity has established that by 2020, forestr y plantations should be managed in ways to ensure the conser vation of biodiversity.This goal aims to r educe pr essur es upon biodiversity, minimizing the ef fects of forestr y practices upon it (CBD2010).To achieve it, forestr y plantations should be a suitable surrogate habitat for native species, and the presence of a well developed understor y is a key factor.Therefore, forestr y plantations holding such undergrowth could provide both goods and services while contributing to conserve a suite of native species (Simonetti et al. 2012).To foster wildlife-friendly practices in other wise "productive lands" (Fisher et al. 2008), an assessment of the understory development is required to evaluate the potential of forestr y stands as surrogate habitat for fauna as a necessary step to integrate these areas devoted to forestry production into wildlife conservation programs as well.The use of LiDAR data in this assessment is feasible, albeit a conser vative approach.The height from DCM had the strongest correlation to understor y heights, but such correlation decreases while the actual height of the dominant canopy, i.e. plantations, increases.This fact could be accounted for by the number of returns available in each vegetation layer, which gets lower when upper canopy layer gets higher.This effect is stronger as plantation density increases (Lefsky et al. 2002).
The predicted understor y height map provided reliable information for ecological applications at a ver y high spatial resolution and it is comparable to similar studies (Turner et al. 2003, Goetz et al. 2007, Vierling et al. 2008, Bergen et al. 2009).Hotspots of higher values of understor y height can be detected and its continuity is captured in the predicted map (Fig. 2 top-right) as can be seen when draw against observed understory height map (Fig. 2 top-left).It is important to consider that estimation were done at one square meter cells and, although simple linear modeling showed low r-square values, the residuals are controlled within the range of 0.24 and 0.41 meters, allowing the use of this type of understor y maps in further studies.However, the estimates exhibit some bias between 0.2 and 0.5 m for high understor y heights when it has values near or above 2.0 m, tending to underestimate their values.This ef fect can be related to the original point cloud acquisition and the vegetation architecture (Lefsky et al. 2002).The gaps in the dominant canopy allow for some of the LiDAR pulses to reach the lower vegetation which has also a low probability to return from the top of it, generating lower heights.As these points can mainly be gather in the gaps, no high angles can be used as they the pulses need to be as vertical as possible to actually receive any returns, which also reduce the possibilities of top understory points to be form in the fi nal cloud.In any event, we can hypothesize that the point cloud density, from underneath the main canopy layer in a given small area, should be related to the understor y coverage and density, and we can only have good samples of this str ucture in the gaps.On another hand, the negative relation between canopy density and understory development has been well documented for managed native forest and plantations (Castedo-Dorado et al. 2012, Shatford et al. 2009, Ogden et al. 1997).This relation states that the lower density in trees the bigger the understory development, which means higher heights.Putting ever ything together, it would be possible to estimate understory density and coverage by using the density of return points in the point cloud or, at least, in DCMs.
One of the advantages of our approach relies on the fact that it requires a minimum of fi eld information, reducing associated time and operational costs.Furthermore, it may be applicable to any type of forests, either native or plantations.This means that is possible to built understor y maps at the landscape level in order to feed the large scales required to properly plan spatially explicit conser vation activities (e.g., Lindenmayer et al. 2006).Another advantage is availability of large LiDAR datasets, increasingly more available for such applications, especially in the central coastal range (e.g., Cartus et al. 2012).In this region, land use change and fragmentations processes are very dynamic (Echeverría et al. 2006) and such maps can play an important role to in biodiversity conser vation planning, potentially incorporation of plantations to ensure both habitat provisioning and landscapelevel connectivity for the fauna.The simple model advanced here will then allow the rapid and conservative estimate of how much of the 1.5 million ha currently allocated to Pinus could be incorporated into wildlife conser vation in Chile as well as abroad, contributing to render the forest industry a more sustainable one.

Fig. 1 :
Fig.1: Study area in Maule Region, South-central Chile.Central positions of clusters are overlaid.The right fi gure represents the inside confi guration of each cluster, and depicts the relative position of sampling subplots.

Fig. 2 :
Fig. 2: Interpolated rasters of both observed (top-left) and predicted (top-right) understory average heights.Bottom left: Map of standard residues of overall linear model.Bottom right: map of slope classes.

Fig. 3 :
Fig. 3: Left: Histogram of residuals from the model that uses all sampling units and all ten predictors.Right: Scatter plot of observed understory heights (H) against predicted ones (fi tted).

TABLE 1
Flight and system variables of the fl ight over Pantanillos Forest Research Station, MauleRegion, South-central Chile.

TABLE 2
Description of variables used in interpolation of heights using LiDAR data available at each stratum of height classes.All variograms were calculated using a lag of 50 m.Descripción de variables empleadas en la interpolación de alturas usando datos LiDAR disponibles de cada clase de altura.Todos los variogramas fueron estimados usando un retardo de 50 m.