Faculty of Science and Forestry
PREDICTION OF SPECIES-SPECIFIC VOLUMES USING SINGLE- SENSOR MULTISPECTRAL AIRBORNE LASER SCANNING AND
NONPARAMETRIC kMSN METHOD Mengqi Li
MASTER’S THESIS
TransAtlantic Forestry Master
JOENSUU 2018
Li, Mengqi 2018. Prediction of species-specific volumes using single-sensor multispectral airborne laser scanning and nonparametric kMSN method. University of Eastern Finland, Faculty of Science and Forestry, School of Forest Sciences.
Master’s thesis in Forest Science, 30 p.
ABSTRACT
Airborne laser scanning (ALS) has been widely used in forest inventory. So far, ALS systems usually comprised of one wavelength (1064 nm or 1550 nm), known as monospectral ALS system.
Teledyne Optech Inc. (Toronto, ON, Canada) introduced the first three-wavelength system, the Titan multispectral airborne lidar, which was equipped with 1550 nm, 1064 nm, and 532 nm lasers.
The objective of this study was to assess the accuracy gains provided by the multispectral ALS data compared to monospectral ALS data in predicting species-specific volumes and classification of dominant species at plot-level. The study area was located in eastern Finland which covered about 428 km2. The species-specific volumes were imputed with kMSN method and area-based approach (ABA) using variables calculated from 1) each wavelength separately, 2) merged wavelength and 3) multispectral model containing all variables calculated from single wavelength and merged wavelength. The predicted species-specific volumes were used for classifying dominant species at plot-level. The best accuracy of species-specific volumes was achieved using multispectral model, resulting in relative root mean square errors (RMSEs) of 24.82%, 63.57%, 53.15%, and 102.42% for total volume, pine, spruce, and deciduous volumes respectively. The results on total volumes among five datasets were close. The accuracies on classifying dominant species were similar between multispectral model and merged wavelength, which were higher than using single wavelength. The improvements from the multispectral model might be due to the multispectral properties of ALS data or due to the higher merged ALS point density. To conclude, I suggest that the multispectral ALS system is further tested and assessed on its potential on providing species-specific volumes using ABA, especially on the effects of increased point density from single wavelength ALS data.
Keywords: species-specific volume, multispectral lidar, airborne laser scanning, kMSN, area- based approach, forestry
1 INTRODUCTION
There has been an increasing demand for accurate and timely information on sustainable forest management around the world. In practices, forest management decisions are usually based on multiple objectives, including timber production, wildlife habitat, biodiversity and carbon balance, which requires information to support decision-making (Temesgen et al. 2007). In Nordic countries, conventional forest inventories at various geographical scales have been enhanced by using remotely sensed data including light detection and ranging (lidar) data, also known as airborne laser scanning (ALS) data, and stereo aerial photography (Anttila 2002, Breidenbach et al. 2010a).
The successful use of ALS data has been demonstrated for a variety of forestry applications. The accurate three-dimensional measurements obtained by ALS data has made ALS technology become one of the most valuable remote sensing methods for providing forest information. By using ALS data, forest attributes such as volume, height, diameter at breast height (DBH), crown area, and stem density can be measured or estimated with high accuracies (i.e. Næsset 2002). As for the use of ALS data in forest inventory, many studies have indicated high accuracy and good usability from an operational perspective (Næsset 1997, 2004, Maltamo et al. 2014).
Tree species information is another important arribute in forest management and planning. Hence, it is necessary to distinguish tree species and obtain information on species-specific attributes (i.e.
volume). However, usually ALS data alone is not suitable for separating species and predicting species-specific attributes (Packalén & Maltamo, 2006, 2007). Most commonly, digital aerial photographs are combined with ALS data (e.g. Packalen & Maltamo 2007). Studies have demonstrated that features of aerial photographs such as tone and texture can be used to discriminate between certain tree species (Franklin, 2001). The near-infrared proportion of the wavelength spectrum (tone) is proved to be especially useful for discriminating between deciduous and coniferous tree species (Kalensky & Wilson, 1975).
ALS data is usually analyzed using two distinguished methods. First one is the individual tree detection (ITD) approach. In this approach, individual tree crowns are segmented either from raster images that derived from canopy height model (CHM) using ALS point clouds (Leckie et al. 2003,
Popescu et al. 2004) or from ALS raw point clouds (Li et al. 2012, Reitberger et al. 2009). In addition, species can be classified at tree-level using shape, structure, and intensity features of tree crowns (Zhang et al. 2015). Tree properties, such as height and biomass, can subsequently be estimated using the segment properties (i.e. crown size) and ALS-derived height metrics as predictor variables (Edson & Ming 2011, Kwak et al. 2007).
The other approach for predicting forest attributes using ALS data is based on height and density variables calculated from ALS data over a certain area, known as the area-based approach (ABA) (Næsset 2002). The ABA was tested experimentally in Finland since 2004 (Maltamo et al. 2006).
Research has shown that ABA can obtain higher accuracy for total stand volume when compared to the previous field-based inventory method or other remote-sensing based inventory method (Uuttera et al. 2006). Later in Finland ABA was modified by Packalén and Maltamo (2006, 2007) to estimate species-specific stand attributes by combining ALS data with aerial imagery, which fulfilled the requirements of stand level management inventory in Finland.
Depending on the properties of ALS data, the purposes of the final products, and stand complexity, both ITD and ABA could be very useful for estimating forest stand attributes. ITD is more intuitive than ABA because the response variables refer to the single tree, which is the smallest unit in forest management. ITD approach is suitable for developing species-specific models, particularly for mixed stands (Kandare et al 2017). In addition, the tree coordinates provided by ITD may be useful in single tree harvesting operations and in growth predictions with distance-dependent growth models (Edson & Wing 2011). However, compared to ABA, ITD usually requires higher density of ALS point clouds data which may lead to higher operational costs. Besides, due to the segmentation errors from ITD, the estimate of the response variable is likely to be biased if aggregated to a larger geographical unit such as a forest stand (Maltamo et al. 2004).
Traditionally, ALS data is monospectral data which is captured using single beam, and the most commonly used wavelength is near-infrared. In December 2014, Teledyne Optech Inc. (Toronto, ON, Canada) introduced the Titan multispectral airborne LiDAR, which was the first commercial three-wavelength system that equipped with scanners emitting three laser beams: 1550nm (middle- infrared), 1064nm (near-infrared) and 532nm (green). The system is a discrete return ALS (up to four discrete XYZ returns per pulse), but each channel can be equipped with a wavelength digitizer to generate full-waveform data. Studies done by Fenandez-Diaz et al. (2016) indicated that its
multispectral capability enabled new applications, including the production of false color active imagery derived from the lidar return intensities and the automated classification of target and land covers. In addition, three wavelengths and three-look angle design provided redundancy and diversity, which were beneficial technically and financially (Fenandez-Diaz et al. 2016). In forestry applications, the commercial multispectral lidar offers the potential to separate woody and leafy canopy elements from non-foliage by utilizing the differences on both structural and reflectance properties (Morsdorf et al. 2009). Furthermore, it has been proposed that multispectral ALS could provide more information on stand vertical structure and species separation than monospectral ALS data (Budei et al. 2017). However, there is one limitation that this scanner has to be operated at a lower flight height than the one just using near-infrared channels because green channel has poorer reflectance from vegetation.
There are several advantages of using ALS multispectral intensities over spectral signatures derived from passive sensors. The measurements of ALS intensity have the advantage of being independent from external illumination conditions. Compared with aerial imaging cases, there would be less variations in intensities caused by variables shadowing (Woodhouse et al. 2011). In addition, tree and ground signals cannot be separated in passive sensors’ measurements of reflectance, whereas the ALS data allows the discrimination by imposing a height threshold.
Nearest neighbour (NN) imputation approaches have been developed and widely used in forestry applications to estimate forest characteristics (LeMay & Temesgan 2005). NN imputation methods are nonparametric in that they do not rely on any underlying probability distribution for estimation (Everitt 1998), and they provide an alternative to traditional parametric regression models (Altman 1992). In NN imputation, the only assumption is that the X-variables have a strong relationship to the Y-variables, and therefore it can be used to impute missing Y-variables (Eskelson et al. 2009).
Sometimes, NN imputation can use X-variables even without a complete knowledge of the complicated relationships between X- and Y-variables (Fehrmann et al. 2008). In NN, the observations with both response (Y) and predictor (X) variables available are references, whereas the observations with only predictor variables available are referred as targets. To improve the prediction accuracy, several close neighbours can be considered instead of using only the closest neighbour. The prediction is then calculated by the weighted mean of the k nearest neighbours (kNN).
The prediction is based on imputing response values of a target to a reference. Usually, the distance measure in NN approaches refer to a thematic distance with respect to observed predictor variables, rather than the geographical distances between units (Everitt 1998). The distance reflects the statistical relationships between the response and predictor variables in the reference dataset (i.e.
the distance is described by the NN model). There are several mathematical methods to determine the proximity between targets and references based on the predictor variables. Two most popular ones are random forests (RF) (Breimann 2001) and most similar neighbour inference (MSN) (Moeur & Stage 1995). In RF, an ensemble of regression trees is grown and the distance between neighbours depends on how often they share the same final node. Examples of using RF include Hudak et al. (2008) who used RF to estimate total basal area and tree density. In MSN, the distance is derived using canonical correlation analysis (Moeur & Stage 1995). Several studies used MSN to predict species-specific volume (Breidenbach et al. 2010b, Kandare et al. 2017, Packalén &
Maltamo 2007).
The first studies using the kNN approach in a remote sensing-based forest inventory context were presented by Kilkki and Paivinen (1987) and Moeur (1987). Since then, kNN methods using remote sensing data have been widely used for forest inventory applications, including the Finnish national forest inventory (Tomppo 1991). Interests have been focused on the prediction of single features such as the basal area, or multivariate attributes, which means predicting several stand characteristics simultaneously (i.e. volume for different species) (Moeur & Stage 1995).
The objective of this study is to evaluate the magnitude of the accuracy gains in the prediction of species-specific volumes provided by the three wavelengths ALS compared to single wavelength ALS. Moreover, the ability of using ALS data alone to classify dominant species at plot-level was examined. So far, the employments of multispectral ALS technologies in forestry applications are quite new. Most of studies that tested the multispectral ALS system have been focused on land cover classification (Hopkinson et al. 2016) and species identification at individual tree level (Budei et al. 2017, Yu et al. 2017), whereas the studies on predicting species-specific volumes using ABA are still missing. To the best of my knowledge, this is among the first attempts to use multispectral ALS data to predict species-specific volumes using the ABA and kMSN method. This study focused on total stand volume, pine, spruce, and deciduous volumes. The general hypothesis is that the combination of ALS data produced in three wavelengths will improve the accuracy on estimating species-specific volumes.
2 MATERIALS AND METHODS
2.1 Study area and field data
The experimental forest area was located in the municipality of Liperi (62°32ʹN 29°23ʹE, UTM zone 35V), which is in the province of Eastern Finland and is part of the North Karelia region (Fig.
1). About 75% of the area was owned by private citizens, and the remainder by government and other forest companies. The size of the study area was about 428 km2. It was a typical Finnish managed boreal forest area, dominated by coniferous species scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies (L.) Karst.). Deciduous trees included downy birch (Betula pubescens Ehrh.) and silver birch (Betula pendula Roth.). The rotation age in Finland is about 80 years, and most of the plantation stands were even-aged, mixed with other tree species.
Figure 1. Location of the study area and placement of field plots.
The field measurements were designed by the University of Eastern Finland (UEF). Plots were placed using systematic cluster sampling (Figure 1). The clusters were placed systematically in the study area and they were 1,200 meters apart from each other. Within each cluster, four field plots were measured, and they were 300 meters apart. A total of 543 sample plots were measured in summer 2016, about half of which were measured by students at UEF and the rest by a local forest company. No stand stratification was used, however the plots placed on seedling stands (i.e. height less than 7 meters) and other open areas were not measured. Among these 543 field plots, there were 391 plots with a size of 254 m2, i.e. they were circular plots with a radius of 9 meters. If the stem density was less than 800 stems per hectare, 12.6 m radius plot was measured with an area of 500 m2.
Scots pine was the main tree species on 52% of the plots and Norway spruce on 37% of the plots.
Deciduous trees, mainly downy birch and silver birch, were usually in the minority of the tree stock. 82% of the plots were considered as one canopy layer and the rest of the plots were multiple layered. The data included both naturally regenerated and plantation forests. The stand development classes consisted of young (24%), mature (52%) and old-growth (24%), while the proportions of the site fertility classes were very rich 2%, rich 15%, medium 48%, rather poor 15%, and poor 1%.
Tree height, DBH, tree class and tree species were measured and recorded for each tree with a DBH greater than 5 cm. Basal area (G) was calculated in each plot and converted to m2 ha-1. The height (hgm) and diameter (dgm) of the basal area median tree for each species in each storey class were calculated in each plot according to the measured tree list. The number of trees was counted in each plot and then scaled up to hectare level to represent stand density (N). Tree volume was calculated using the species-specific models presented by Laasasenaho (1982) using DBH and tree height as independent variables. Finally, total tree volume (V) was scaled up to volume per hectare and summarized for each plot, all by tree species (i.e. pine, spruce and deciduous trees). Some of the key stand characteristics are presented in Table 1. The minimum and maximum total volume at plot-level are 16.6 and 1012.6 m3 ha-1, respectively.
Table 1. Means and standard deviations (in parentheses) of stand characteristics by tree species at plot level.
V, m3 ha-1 G, m2 ha-1 N, stem ha-1 dgm, cm hgm, m
Pine 74.5 (87.6) 8.4 (9.0) 332 (425) 14.2 (1.2) 11.7 (9.7)
Spruce 86.8 (106.2) 9.8 (10.2) 526 (507) 15.8 (9.8) 13.4 (7.8)
Deciduous 26.9 (42.7) 3.4 (4.9) 261(395) 11.9 (9.6) 12.0 (8.4)
Total 188.2 (105.9) 21.6 (8.7) 1118 (640) - -
2.2 Airborne laser scanning data
The ALS data for the Liperi forest area were collected on 2nd, 3rd, and 10th July 2016, using the Optech Titan multispectral laser scanning system mounted on a fixed wing airplane and operated at an altitude of 900 m above ground level (a.g.l). The Titan system was comprised of lasers firing at three different angles (Table 2). There were three laser beams used: 532 nm (green), 1064 nm (near infrared) and 1550 nm (short-wave infrared). The flight speed was 77 m s-1 and a pulse frequency was 250 kHz. The resulted strip width on average was 655 m. Altogether 82 laser strips were acquired with an overlap of 55% with both first and last return data were recorded. The area covered by the ALS data was approximately 428 km2.
Table 2. Titan system's laser characteristics.
Channel Wavelength Divergence (1/e) Forward tilt Pulse densitya Pulse width
C1 1550 nm 0.35 mrad 3.5° 4.8 pulse/m2 3.0-3.5 ns
C2 1064 nm 0.35 mrad 0.0° 4.8 pulse/m2 3.0-3.5 ns
C3 532 nm 0.7 mrad 7.0° 3.7 pulse/m2 2.5-3.0 ns
a Dr Petteri Packalén, University of Eastern Finland, personal communication.
2.3 Pre-processing of ALS data
Recorded intensity is the amount of energy reflected back to the laser sensor, which is affected by several factors, including target surface characteristics (reflectance, wetness and roughness), environmental effects (atmospheric transmittance, moisture), and acquisition parameters and instruments (Korpela et al., 2010). Therefore, intensity calibration is an important step to compensate the impact of these factors and achieve better classification accuracy. In this study, the ALS intensity values were calibrated for different scan ranges by using the equation described by Korpela et al. (2010):
In = (R/Rref)αIraw,
where In is the calibrated intensity, Iraw is the raw intensity, R is the range from the sensor to targets, Rref is the reference range or average flying height (in this study Rref = 830 m), and α is the optimized power. The optimization was done so that the plot level mean intensities would be as similar as possible when the plots were observed from several scan lines and have different distances. In this study, the correction parameters α were optimized separately for each wavelength. For 1550 nm, 1064 nm, and 532 nm, the α values were 1.1, 0.4, and 0.0, respectively (P. Packalén, personal communication, February 16th, 2017).
The digital terrain model (DTM) was constructed separately for each channel using the ALS points only from that channel. However, since the green light (532 nm) would be absorbed strongly by vegetation, the 532 nm ALS data was normalized using DTM extracted from 1064 nm ALS data.
The first step was to separate the ground laser points from non-ground hits, following the methods explained by Axelsson (2000). Then, a raster DTM with a 1-m pixel size was created from the ground returns by taking the average height of the ground returns within a grid. Canopy height model (CHM) was calculated by subtracting the DTM at the corresponding point for each wavelength. The ALS hits that were less than 2 meters were omitted from analysis to reduce the random variation in heights and the effects of understory vegetations. Two datasets were prepared for each wavelength: one includes only first and single returns and the second includes only last and single returns. These two datasets were spatially registered to the field plots using software package FUSION. Pulses that hit outside these field plots were excluded from further analysis.
2.4 ALS metrics
The ALS predictor variables used in the MSN imputation were derived from the ALS height distributions. Calculations was undertaken using built-in function cloudmetrics provided by FUSION software package. Point density was calculated for each wavelength separately (Table 3).
All the ALS metrics were calculated twice using only first and single returns and only last and single returns (Table 4). The calculations were done for each wavelength and for the merged wavelength which combined ALS point clouds from three wavelengths together. In addition, a multispectral dataset was created which included all ALS-derived variables from three monospectral ALS data and the merged ALS data. This resulted in five ALS-derived variable datasets, named as Wavelength 1, Wavelength 2, Wavelength 3, Merged Wavelength, and Multispectral Model. For the first four datasets, the variables were the same except that they were
calculated from different wavelength or merged wavelength, resulting in 62 variables for each dataset. The last dataset included all variables from individual dataset, plus the intensity ratio variables, resulting in 252 variables in total. The details of calculated variables are summarized in Table 4.
Table 3. ALS point density statistics calculated from first returns data for each wavelength and the merged data.
ALS point density statistics (points m-2)
25th percentile Median 75th percentile Mean
1550nm 9.8 11.2 18.0 14.8
1064 nm 10.0 11.2 17.6 14.7
532 nm 5.9 7.7 10.4 7.7
Merged wavelength 26.5 30.2 45.8 38.3
Table 4. Stand features derived from normalized point data and spectral information.
ALS-derived variables Definition Point cloud variables
Havg Arithmetic mean of normalized height of all points above 2 m threshold Hstd Standard deviation of normalized height of all points above 2 m threshold
Hmax Maximum of the normalized heights of all points
cc_2 Percentage of points above 2 meters indicating crown cover
cc_avg Percentage of points above the average height indicating crown cover HP10 to HP99 10% to 90% percentiles of normalized height of all points with a 10%
increment Single-channels intensity variables
Int_avg Mean of intensity
Int_std Standard deviation of intensity
Int_max Maximum of intensity
Int_sk Skewness of intensity
Int_kut Kurtosis of intensity
Int10 to Int99 Percentiles of intensity from 10% to 90% with 10% increments Multi-channel intensity variablesa
Int_b2b1 The ratio of average intensity between channel 2 and channel 1 Int_b3b1 The ratio of average intensity between channel 3 and channel 1 Int_b3b2 The ratio of average intensity between channel 3 and channel 2 Ratio_b2b1 The normalized difference (Int_b2 – Int_b1)/( Int_b2 + Int_b1)
a These variables were only available in the Multispectral Model dataset.
2.4 Estimation of stand characteristics 2.4.1 kMSN imputation
Stand attributes of interest at plot level included total volume (m3 ha-1) and volumes (m3 ha-1) by main tree species: pine, spruce, and deciduous. The stand attributes were predicted simultaneously by nonparametric k-nearest neighbor method utilizing k most similar neighbor (kMSN) imputation.
The kMSN method is a modification of the original MSN method proposed by Moeur and Stage (1995). It uses canonical correlation analysis to produce a weighting matrix for selecting the k most similar neighbours from the reference data. Canonical correlations are implemented to find the linear transformations Ul and Vm for the set of dependent variables (Y) and predictor variables (X) which maximize the correlation between them:
Ul = αlY, and Vm = γmX,
where αl represents the canonical coefficients of the dependent variables, γm represents the canonical coefficients of the predictor variables, and l and m represent the number of Y and X variables, respectively. Ul and Vm are ordered in such a way that the canonical correlation is largest for l or m = 1, second for l or m = 2, etc. Thus, the predictive power is concentrated in the first few canonical components.
The MSN distance metric derived from the canonical correlation analysis is as follows (Moeur &
Stage 1995):
where Xu is the vector of the predictor variables from the target observation, Xj is the vector of the predictor variables from the reference observation, Γ is the matrix of canonical coefficients of the predictor variables and Λ is the diagonal matrix of squared canonical correlations. The kMSN method used here is the same as MSN except that estimates for the target observation are calculated as weighted averages of the k nearest observations and the weighting is based on the inverse of the MSN distance. The weighted Wuj of a reference plot u for the target plot j was calculated as follows where k is the number of nearest observations:
.
Nearest neighbour search and imputation was done by using yaImpute R package developed by Crookston and Finley (2008). Two key functions used in yaImpute package are impute and yai.
They were used in searching for the best combinations of predictor variables and imputation processes.
2.4.2 Predictor variable selection
The selection of predictor variables has been found to be laborious and time-consuming, and therefore should be automated (Maltamo et al. 2006, Packalén & Maltamo 2006). Selection of predictor variables is necessary for nonparametric approaches as it can reduce the computation time and solve the collinearity. In this study, variable selection for the modelling was done using an algorithm that implemented optimization based on simulated annealing in R. The R package GENSA was used in the algorithm, to minimize RMSE% for the combination of four dependent variables. The weights for four response variables were equal. During preliminary tests, different numbers of predictor variables (x) and neighbours (k) were tested in order to get the optimal combinations of x and k.
For all five datasets, different numbers of predictor variables (x) were chosen (i.e., x = 15, 20, 25) using the automatic selection algorithm. Based on the preliminary tests and literatures, the number of predictor variables was set to 20 for all five datasets. However, the slight differences on x (i.e., 19 vs. 21) only affected the imputation accuracy slightly with respect to the response’s RMSE%
values. The transformation of the predictor variables was not applicable in this study as the main purpose for this study was to compare the suitability of using multispectral ALS with monospectral ALS data. Since the transformation only improved the model accuracy marginally and the number of candidate variables was large, they were not included here.
During the preliminary tests and the literature reviews, it was noted that for all imputations, as k (number of neighbours) increased, the RMSE% values decreased dramatically at the beginning (i.e., k = 1 to 5), levelled out to a value of approximately 7, and started to increase again for large k values (i.e., k > 30). In this study, k = 10 was chosen for all imputations and for further analysis.
Once the number of predictor variables and the number of neighbours were determined, the predictor variables (Table 4) were selected separately for each wavelength using the automatic selection algorithm, so that the lowest average RMSE% of response variables can be achieved.
2.5 Dominant species classification
In addition to the species-specific volumes, each plot was also classified according to its main tree species to test the suitability of dominant species classification using multispectral ALS data. Once the species-specific volume was predicted for each plot, the species with the highest volume was identified and treated as the dominant tree species for that plot. The similar procedure was done for the reference data using the original field measured volumes for each species. The results of species classification at plot level were then compared between the predicted and reference records. The same process was done using predicted volumes from each ALS dataset.
2.6 Accuracy assessment
Results were validated by means of leave-one-out cross-validation. Reliability characteristics were calculated for each predicted variable by tree species and for the total volume. The accuracy of the estimates is expressed in terms of the root mean square error (RMSE):
RMSE = %"&'!"#!" $ ( ,
where n is the number of sample plots, yi is the observed value for plot i and 𝑦* is the predicted value for plot i. Moreover, the relative RMSEs (RMSE%) were calculated by dividing the absolute RMSE values by the true mean for the characteristic concerned. Additionally, the squared correlation (r2) was computed as the Pearson’s correlation coefficient of observed and predicted values. The biases were calculated as follows:
bias = %"&'(!"#!") ( .
To assess whether the bias was statistically significantly different from zero, the following test statistic was used:
tbias = // (-*./ ,
where s is the standard deviation of the differences between ground truth observations and their corresponding estimates.
The accuracy of dominant species classification was evaluated by comparing the classified dominant species with the reference stand calculated using data recorded in the field. The result of the comparison can be represented by an error matrix. Four widely used measures, producer’s accuracy, user’s accuracy, overall accuracy (OA), and Kappa coefficient, were computed to evaluate the performance of the classification.
3 RESULTS
3.1 Variable and parameter selection
Table 5 shows that for the monospectral ALS datasets, the number of height-related predictor variables was generally the same as the intensity-related features. For example, for wavelength 1, wavelength 2, and wavelength 3, there were 8, 10, and 10 intensity-related variables, respectively.
For the multispectral model, only 6 height-related features were selected, and 14 features were intensity-related.
Table 5. Predictor variables from ALS derived metrics selected for each dataset.
Dataset Predictor selecteda 1550 nm
(C1)
cc_f2, cc_favg, h_lavg, h_fstd, h_f30, h_l50, h_f60, h_l70, h_f90, Int_lmax, Int_fsk, Int_lavg, Int_fkurtosis, Int_f10, Int_l20, Int_l40, Int_f50, Int_l70, Int_f80, Int_f90, canopy_ratio_l 1064 nm
(C2)
cc_f2, cc_l2, cc_lavg, cc_favg, h_fkurtosis, h_fstd, h_lstd, h_l20, h_l30, h_f80, h_f90, Int_lavg, Int_fmax, Int_lmax, Int_fsk, Int_lsk, Int_lkurtosis, Int_f20, Int_f50, canopy_ratio_l
532 nm (C3)
cc_lavg, h_lavg, h_fstd, h_lsk, h_fkurtosis, h_f20, h_f50, h_l80, Int_favg, Int_fmax, Int_lkurtosis, Int_f10, Int_l10, Int_f30, Int_l50, Int_f70, Int_l70, Int_l80, canopy_ratio_f, canopy_ratio_l
Merged data cc_l2, cc_lavg, h_lmax, h_l10, h_l70, h_f80, h_f90, Int_lavg, Int_fstd, Int_lstd, Int_lkurtosis, Int_lsk, Int_lmax, Int_f10, Int_l10, Int_l30, Int_f60, Int_f70, Int_f90, canopy_ratio_l
Multispectral modelb
cc_favg1, h_lstd2, h_lmax1, f_h402, h_l90123, Int_fmax123, Int_fmax1, Int_fsk2, Int_fstd1, Int_fstd2, Int_f10123, Int_l402, Int_f60123, Int_f603, Int_l603, Int_f701, Int_f80123, Int_l902, canopy_ratio_l3, Int_b2b1
a f and l indicated if the predictor was calculated using first and single returns or last and single returns.
b The subscript indicated which wavelength that the variable came from, where 1, 2, 3 and 123 denoted wavelength 1 (1550 nm), wavelength 2 (1064 nm), wavelength 3 (532 nm) and the merged data, respectively.
3.2 Accuracy of stand attributes
The species-specific volumes and the total volume were imputed using selected predictor variables (Table 5) and k = 10. Figure 2 to 6 show the observed versus predicted volume for each species and the total volume. Table 6 summarizes the RMSE, RMSE%, r2, and bias for each predicted volume. Using t-test and alpha level of 0.05, none of the bias was significant. In general, the predicted total volume had the highest accuracy, followed by spruce volume, pine volume, and the deciduous volume had the lowest accuracy. In terms of the absolute errors (RMSE), total volume, pine and spruce volume had a similar level, whereas deciduous species had the lowest error level because of its minority among other species.
Figure 2. Observed versus predicted total volume and species-specific volumes at the plot level using features derived from Wavelength 1.
Figure 3. Observed versus predicted total volume and species-specific volumes at the plot level using features derived from Wavelength 2.
Figure 4. Observed versus predicted total volume and species-specific volumes at the plot level using features derived from Wavelength 3.
Figure 5. Observed versus predicted total volume and species-specific volumes at the plot level using features derived from Merged Wavelength.
Figure 6. Observed versus predicted total volume and species-specific volumes at the plot level using features from Multispectral Model.
Table 6. Accuracy of the estimated species-specific volume and total volume using different ALS data.
RMSE (m3/ha) RMSE% r2 BIAS Total Volume
Wavelength 1 47.01 24.98 0.81 2.19
Wavelength 2 50.21 26.68 0.78 2.94
Wavelength 3 49.14 26.11 0.78 1.82
Merged Wavelength 48.13 25.57 0.80 2.97
Multispectral Model 46.71 24.82 0.81 1.05
Pine Volume
Wavelength 1 56.99 76.53 0.58 2.59
Wavelength 2 53.26 71.52 0.63 1.78
Wavelength 3 64.09 86.07 0.46 0.92
Merged Wavelength 55.13 74.03 0.61 2.90
Multispectral Model 47.34 63.57 0.71 1.94
Spruce Volume
Wavelength 1 54.50 62.80 0.74 -0.98
Wavelength 2 50.84 58.58 0.77 0.77
Wavelength 3 64.75 74.62 0.63 1.01
Merged Wavelength 48.75 56.17 0.79 -0.07
Multispectral Model 46.14 53.15 0.81 -2.13
Deciduous Volume
Wavelength 1 33.23 123.34 0.39 0.58
Wavelength 2 31.96 118.64 0.43 0.39
Wavelength 3 38.38 142.47 0.20 -0.12
Merged Wavelength 30.30 112.48 0.50 0.14
Multispectral Model 27.59 102.42 0.59 1.23
None of the biases was significant at the 95% confidence level (tbias<1.96).
Among three monospectral ALS data, Wavelength 2 had the lowest averaged RMSE% for species- specific volumes and Wavelength 3 had the worst performance on predicting species-specific volumes. For all five imputations, the RMSE% varied from 24.8% to 26.7% for the total volume, 63.6% to 86.1% for the pine volume, 53.2% to 74.6% for the spruce volume, and 102.4% to 142.5%
for the deciduous volume. It was noticed that using the multispectral dataset generated the lowest RMSE% for total volume (24.8%), pine (63.6%), spruce (53.2%), and deciduous volume (102.4%).
The merged wavelength also did a better job on predicting spruce and deciduous volume than monospectral ALS data.
3.3 Dominant species classification with different wavelengths
The confusion matrix of dominant species classification at plot-level are presented in Table 7-11 based on the kMSN imputations from different wavelengths, i.e. 1,550 nm (Table 7), 1,064 nm (Table 8), 532 nm (Table 9), merged data (Table 10), and multispectral model (Table 11). The highest level of accuracy (84.35%) was obtained with the multispectral model. Using features from 532 nm produced the lowest overall accuracy of 75.32%. It was also noted that the differences on
OA between 1550 nm (80.11%) and 1064 nm (81.77%) were marginal, and it was the same for merged data (83.79%) and multispectral model (84.35%). When using features from three wavelengths together, either using metrics from merged point clouds or all features combined, the performance on stand classification was better than using monospectral ALS data.
Table 7. Confusion matrix and accuracy evaluation of stand classification using ALS metrics from Wavelength 1.
Predicted Producer (%)
Pine Spruce Deciduous
Reference Pine 212 33 6 84.46
Spruce 16 215 4 91.49
Deciduous 18 31 8 14.04
User (%) 86.18 77.06 44.44 OA = 80.11%, Kappa = 0.65
Table 8. Confusion matrix and accuracy evaluation of stand classification using ALS metrics from Wavelength 2.
Predicted Producer (%)
Pine Spruce Deciduous
Reference Pine 225 21 5 89.64
Spruce 17 210 8 89.36
Deciduous 21 27 9 15.79
User (%) 85.56 81.40 40.91 OA = 81.77%, Kappa = 0.68
Table 9. Confusion matrix and accuracy evaluation of stand classification using ALS metrics from Wavelength 3.
Predicted Producer (%)
Pine Spruce Deciduous
Reference Pine 210 39 2 83.67
Spruce 31 197 7 83.83
Deciduous 21 34 2 3.51
User (%) 80.15 72.96 18.18 OA = 75.32%, Kappa = 0.56
Table 10. Confusion matrix and accuracy evaluation of stand classification using ALS metrics from Merged Wavelength.
Predicted Producer (%)
Pine Spruce Deciduous
Reference Pine 222 23 6 88.45
Spruce 18 213 4 90.64
Deciduous 15 22 20 35.09
User (%) 87.06 82.56 66.67 OA = 83.79%, Kappa = 0.72
Table 11. Confusion matrix and accuracy evaluation of stand classification using ALS metrics from Multispectral Model.
Predicted Producer (%)
Pine Spruce Deciduous
Reference Pine 223 22 6 88.84
Spruce 14 211 10 89.79
Deciduous 17 16 24 42.11
User (%) 87.80 84.74 60.00 OA = 84.35%, Kappa = 0.73
Classification accuracies also varied among species. The highest accuracies were obtained for pine stand with an 89.64% producer’s accuracy using ALS metrics from 1,064 nm, and for spruce stand
with a 91.49% producer’s accuracy using ALS metrics from 1,550 nm, followed by deciduous species with 42.11% producer’s accuracy using all features combined data (multispectral model).
When comparing among 1,550 nm, 1,064 nm, merged data, and all features combined data, it was noted that the producer’s accuracies for classifying pine and spruce stands were quite similar (i.e.
ranges from 84.46% to 88.84% for pine stands and 89.36% to 91.49% for spruce stands), whereas it was the deciduous stand classification which has the highest variation on producer’s accuracies (ranged from 14.04% to 42.11%).
4 DISCUSSION
This study demonstrates the use of kMSN methods and multispectral ALS to separate dominant species and estimate species-specific volumes and compare the results to that obtained from using monospectral ALS data. Results suggested that the multispectral model, which used ALS metrics from merged point clouds and monospectral ALS data, had the lowest RMSE% for total volume and species-specific volumes. This improvement might due to the multispectral properties of the data. For example, some of the features from merged data were helpful on separating species, and when they were combined with features from monospectral ALS data, they increased the general accuracies for all species as well as the total volume. However, the improved accuracy might be due to the increased ALS point density from merged data instead of the multispectral properties.
In other studies, preliminary tests were done using density-normalized point data to assess the effect of point density on predicting forest biomass using area-based approach, and they showed that the differences on accuracies between full data and normalized data were minor (M. Maltamo
& L. Korhonen, personal communication, April 3rd 2018). Further studies should be carried out using density-normalized monospectral and merged ALS data to reduce any potential effects of point density on the model accuracy.
The dominant species classification has better accuracy on coniferous plots compared to deciduous plots. However, it should be noted that this classification was done using the species with highest volume, and it did not show the species composition if it was a mixed-species plot (the volume proportions of three species groups were relatively close, i.e. 33%, 33%, and 34% for pine, spruce
and deciduous). In that case, the plot should be classified as mixed-species rather than single- species dominant plot, and this would affect the overall classification accuracy of this study.
It can be problematic to compare the accuracies of this study directly to other studies using multispectral ALS data because of the differences in the data used, the number and type of species identified, and the complexity of stand structure. The relative RMSEs on species-specific volumes achieved in this study were similar to the study done by Packalén and Maltamo (2007) who used the combination of ALS data and aerial images and kMSN area-based approach. Both study areas were located in southern Finland where the species composition and stand structure were similar.
In their best models, the RMSE% for total volume, pine, spruce, and deciduous at plot level were 20.51%, 51.55%, 55.72%, and 102.84%, respectively. In this study, the lowest RMSE% for these four response variables were 24.82%, 63.57%, 53.15%, and 102.42%, respectively. One possible reason for the lower RMSE% on total and pine volume in Packalén and Maltamo’s study might be that they included metrics from aerial images in addition to ALS metrics. Moreover, the study area in their study also has less variation in forest structures because the similar forest silviculture practices were applied to the entire area.
Using different plot size could be a concern in area-based approach because it would lead to bias on estimated volumes and affect the evaluation on the performance of monospectral and multispectral ALS data. Larger plot size generally had better accuracy on estimating volumes, although the exact effects on the results (i.e. how much bias) for this study were unknown. The best way to reduce this unknown bias is to use the same plot size throughout the entire study area in future studies.
In this study, deciduous species were grouped together as it was suggested that their separation by remote sensing methods is challenging under Finnish forest conditions, where the forests are dominated by coniferous species and the vast majority of the deciduous trees are birches (Packalén
& Maltamo 2007). In addition, deciduous trees often grow below the dominant coniferous tree layer. The absolute errors caused by deciduous trees are quite small because of their small quantities in Finnish forests. However, in other countries where the species composition is more complex, or the separation of deciduous trees is important, the suitability of using multispectral ALS and ABA is unknown and should be studied for more details.
Intensity information was important in separating species in this study. Not only the ratio between averaged intensity from different wavelengths, but also different intensity percentiles from different wavelengths were used in the multispectral imputation, and the total number of intensity- related variables were 14 out of 20 variables, indicating that intensity variables may be more useful on separating species than height variables. Korpela et al. (2010) classified Scots pine, Norway spruce, and birch by using intensity variables at Hyytiälä in southern Finland and showed that intensity features can contribute to a classification accuracy of 88% among these three species.
Although the quality of intensity values extracted from multispectral ALS data was not tested or compared to those from fused ALS data, other studies showed that the use of single source ALS data has advantages over the use of fused ALS data in terms of data processing (Puttonen et al.
2010).
Multispectral ALS data was the only remotely-sensed data source used in this study. However, it may be worthwhile investigating whether multispectral ALS data should be combined with other sources of remotely-sensed data (i.e. satellite images) in order to improve the accuracies on species- specific estimates. In addition, only volumes were tested in this study, and the performance of multispectral data on predicting other features (i.e. basal area) was unknown. ALS data alone may not very well suited for the estimation of other species-specific forest characteristics at the plot level. Different remotely sensed technologies sense different aspects of forest structure, so the integration of data from different remote sensors may be worthwhile in providing useful and relevant information (Hudak et al. 2002, LeMay et al. 2008, Törmä 2000). For example, features of aerial photographs such as tone and texture can be used to discriminate between certain tree species (Franklin 2001). Landsat imagery is useful for classifying the spatial extent and seasonal phenology of forest stands across a landscape.
Geographical variables may also be considered in further studies. It was found that topographical parameters (i.e. aspect and elevation) (Hudak et al. 2008) and site characteristic data (i.e. site fertility and ecological zone) can be very helpful to derive species specific information, although other studies (Kotivuori et al. 2018) argued that geographical information may be more useful for much larger areas (i.e. national scale). However, such parameters were not included in this study for these reasons: 1) the study area was quite homogenous on these characteristics (i.e. aspect and elevation), and there was not enough variability so that the contributions from these geographical layers were limited. 2) There was only one class (i.e. ecozone) for the whole study area, so this
information was not useful. 3) Based on the preliminary tests, the relationships between some of the geographical information (i.e. site fertility class) and species were not strong. To reduce the model complexity, these variables were not included in the analysis. Besides, these variables are not always available for the study area.
The use of models outside the range of the modeling data may result in unreasonable estimates.
This study only included the analysis of trees with DBH greater than 5 cm and did not include any smaller trees. Thus, the imputation would not be suitable for predicting young stands or seedlings.
The accuracy of kMSN imputation would be dependent on the variability of reference data (Moeur
& Stage 1995), which means that it may not work well for extreme values (i.e. very high or low volume) where there is not enough reference data to search for k nearest neighbours for the imputation. For example, there was one plot with total volume greater than 1,000 m3 ha-1 which had relatively poor estimation, because that was the only plot with high volume and there were no neighbours available. The accuracy of the estimates improves with increasing k to an optimal choice of k. However, it is found that the estimates may not be within the bounds of reality if more than one neighbour is used (LeMay & Temesgen 2005), and the estimation precision for extreme values of Y-variables decreases with an increase in k (McRoberts et al. 2002).
During preliminary tests, it was noticed that different predictor variables were identified during different runs of the automated variable selection procedure, which was also reported by Packalén et al. (2009). This is caused by the fact that many predictor variables are highly correlated, and the same variable selection is consequently not possible based on the data alone (Faraway 2002). The choice of X-variables depends on the information that is available and on the variables related to the Y-variables (LeMay & Temesgen 2005). Increasing the number of X-variables does not guarantee improvement in the estimation results (McRoberts et al. 2002). As the number of X- variables increases, it becomes increasingly difficult to find relevant neighbours (Maltamo et al.
2006). The optimal value for k is a trade-off between the accuracy of the estimates and the variation that is retained in the estimates (McRoberts et al. 2002), which complicates the use of nonparametric methods. The choice of optimal value of k is also affected by the strength of the relationship between the X- and Y-variable, meaning that with weaker relationships would result in larger k values (McRoberts et al. 2002).
The weights for each response variable are assumed to be equal in the variable automatic search algorithm. Another alternative could be giving different weights to each response variable. Usually, the variable selection algorithm is designed to minimize the weighted average of the RMSEs, and it makes sense to give more weight to those response variables which require most improvements, whereas for other dependent variables, the accuracies are almost always at an acceptable level (Packalén & Maltamo 2007). The selection of weights, however, is somewhat subjective. In this study, the main purpose is to compare the accuracies of using different sets of ALS data, so as long as the weight vectors for each response variable are the same, the results are comparable. The optimal choice of k, the distance metric including weights, and X-variables is difficult to determine (LeMay & Temesgen 2005). The best combination depends on the problem and the available data.
In this study, ALS data was collected in summer, which is also known as leaf-on data. As an alternative to leaf-on data, there were some studies showing that leaf-off data may be helpful to discriminate coniferous from deciduous species and improve the accuracy on species-specific estimations, although it may depend on the characteristics of forest stands. (Villikka et al. 2012).
In the future, multispectral ALS data could be collected during leaf-off season to test how this affect the accuracies on the species-specific estimates. However, there are also some issues with leaf-off data according to previous studies with monospectral ALS data. First, leaf-off ALS data tends to underestimate the tree height than leaf-on data, especially in the case of deciduous species (Gaveau & Hill 2003), as under leaf-off conditions, the laser pulses penetrate deeper into deciduous tree crowns than under leaf-on conditions before an echo is detected. Second, the timing of gathering leaf-off data is critical in Finland, as the response is affected by snow cover on the ground. Since the time of snow melting varies from year to year, it may be problematic to gather data during short time window (Liang et al. 2007, Kim et al. 2009).
5 CONCLUSION
In this study, I examined the potential utility of single-sensor multispectral ALS data for estimating species-specific volumes in a conifer-dominant stand in a boreal forest zone. The results suggest that additional information, provided by multispectral laser scanning, may be a valuable source of information for separating species and providing species-specific volumes on pine, spruce, and
deciduous trees, which are the main tree species (species group) found in boreal forest zones.
However, further studies are needed on the effects of ALS point densities from monospectral and multispectral data on the accuracy of estimates. Currently, it is more expensive to acquire multispectral ALS data than aerial images and single-channel ALS data. It is anticipated that the price will drop as technology develops, and the market is growing. Before expanding the operation of using multispectral ALS data as solution for estimating species-specific volumes, it requires more detailed studies and tests on different study areas.
REFERENCE
Altman, N. S. (1992). An introduction to kernel and nearest-neighbor nonparametric regression.
The American Statistician, 46(3), 175-185. 10.1080/00031305.1992.10475879
Anttila, P. (2002). Nonparametric estimation of stand volume using spectral and spatial features of aerial photographs and old inventory data. Canadian Journal of Forest Research, 32(10), 1849-1857. 10.1139/x02-108
Axelsson, P. (2000). DEM generation from laser scanner data using adaptive TIN models. Proc.
XIXth ISPRS Congress, IAPRS, Vol. XXXIII. Amsterdam, The Netherlands.
Breidenbach, J., Nothdurft, A., & Kändler, G. (2010a). Comparison of nearest neighbour approaches for small area estimation of tree species-specific forest inventory attributes in central Europe using airborne laser scanner data. European Journal of Forest Research, 129(5), 833-846. http://dx.doi.org/10.1007/s10342-010-0384-1
Breidenbach, J., Næsset, E., Lien, V., Gobakken, T., & Solberg, S. (2010b). Prediction of species specific forest inventory attributes using a nonparametric semi-individual tree crown approach based on fused airborne laser scanning and multispectral data. Remote Sensing of Environment, 114(4), 911-924. 10.1016/j.rse.2009.12.004
Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5-32. 10.1023/A:1010933404324 Budei, B.C., St-Onge, B., Hopkinson, C., & Audet, F. (2017). Identifying the genus or species of
individual trees using a three-wavelength airborne LiDAR system. Remote Sensing of Environment, 204, 632-647.10.1016/j.rse.2017.09.037
Crookston, N. L., & Finley, A. O. (2008). yaImpute: An R Package for k-NN imputation. Journal of Statistical Software, 23(10), 1-16.
Eskelson, B. N. I., Temesgen, H., Lemay, V., Barrett, T. M., Crookston, N. L., & Hudak, A. T.
(2009). The roles of nearest neighbor methods in imputing missing data in forest inventory and monitoring databases. Scandinavian Journal of Forest Research, 24(3), 235-246.
10.1080/02827580902870490
Edson, C., & Wing, M. (2011). Airborne light detection and ranging (LiDAR) for individual tree stem location, height, and biomass measurements. Remote Sensing, 3(11), 2494-2528.
doi:10.3390/rs3112494
Everitt, B. S. (1998). The Cambridge dictionary of statistics. Cambridge: Cambridge University Press.
Faraway, J.J. (2002). Practical regression and anova using R. Ann Arbor, MI, self-published.
http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf.
Fehrmann, L., Lehtonen, A., Kleinn, C. & Tomppo, E. (2008). Comparison of linear and mixed- effect regression models and a k-nearest neighbour approach for estimation of singletree biomass. Canadian Journal Forest Research, 38, 1-9.
Fernandez-Diaz, J.C., Carter, W.E., Glennie, C., Shrestha, R.L., Pan, Z., Ekhtari, N., … & Sartori, M. (2016). Capability assessment and performance metrics for the titan multispectral mapping Lidar. Remote Sensing, 8(11), 936.10.3390/rs8110936
Franklin, S.E. (2001). Remote sensing for sustainable forest management. Boca Raton: CRC Press LLC.
Gaveau, D., & Hill, R. (2003). Quantifying canopy height underestimation by laser pulse penetration in small-footprint airborne laser scanning data. Canadian Journal of Remote Sensing, 29(5), 650-657. 10.5589/m03-023
Hopkinson, C., Chasmer, L., Gynan, C., Mahoney, C., & Sitar, M. (2016). Multisensor and multispectral LiDAR characterization and classification of a forest environment. Canadian Journal of Remote Sensing, 42(5), 501-520. 10.1080/07038992.2016.1196584
Hudak, A. T., Crookston, N. L., Evans, J. S., Hall, D. E., & Falkowski, M. J. (2008). Nearest neighbor imputation of species-level, plot-scale forest structure attributes from LiDAR data. Remote Sensing of Environment, 112(5), 2232-2245. 10.1016/j.rse.2007.10.009 Immitzer, M., Atzberger, C., & Koukal, T. (2012). Tree species classification with random forest
using very high spatial resolution 8-band WorldView-2 satellite data. Remote Sensing, 4(9), 2661-2693. 10.3390/rs4092661
Kalensky, Z., & Wilson, D. A. (1975). Spectral signatures of forest trees. Proc. 3rd Canadian Symposium on remote sensing, Edmonton, Alberta, 22 -24 September.
Kandare, K., Dalponte, M., Orka, H., Frizzera, L., & Naesset, E. (2017). Prediction of species- specific volume using different inventory approaches by fusing airborne laser scanning and hyperspectral data. Remote Sensing, 9(5), 400. 10.3390/rs9050400
Kilkki, P. & Paivinen, R. (1987). Reference sample plots to compare field measurements and satellite data in forest inventory. In Proceedings of Seminar on remote sensing-aided forest inventory, Hyytiala, Finland. Vol. 19. Department of Forest Mensuration and Management, Research Notes (pp. 209-215). Helsinki: University of Helsinki.
Kim, S., McGaughey, R. J., Andersen, H., & Schreuder, G. (2009). Tree species differentiation using intensity data derived from leaf-on and leaf-off airborne laser scanner data. Remote Sensing of Environment, 113(8), 1575-1586. 10.1016/j.rse.2009.03.017
Korpela, I., Orka, H., Maltamo, M., Tokola, T., & Hyyppa, J. (2010). Tree species classification using airborne LiDAR - effects of stand and tree parameters, downsizing of training set, intensity normalization, and sensor type. Silva Fennica, 44(2), 319-339. 10.14214/sf.156 Kotivuori, E., Maltamo, M., Korhonen, L., & Packalen, P. (2018). Calibration of nationwide
airborne laser scanning based stem volume models. Remote Sensing of Environment, 210, 179-192. doi:10.1016/j.rse.2018.02.069
Kwak, D., Lee, W., Lee, J., Biging, G. S., & Gong, P. (2007). Detection of individual trees and estimation of tree height using LiDAR data. Journal of Forest Research, 12(6), 425-434.
doi:10.1007/s10310-007-0041-9
Laasasenaho, J. (1982). Taper curve and volume function for pine, spruce and birch: Männyn, kuusen ja koivun runkokäyrä- ja tilavuusyhtälöt. Helsinki: Finnish Forest Research Institute.
Leckie, D., Gougeon, F., Hill, D., Quinn, R., Armstrong, L., & Shreenan, R. (2003). Combined high-density lidar and multispectral imagery for individual tree crown analysis. Canadian Journal of Remote Sensing, 29(5), 633-649. 10.5589/m03-024
LeMay, V., & Temesgen, H. (2005). Comparison of nearest neighbor methods for estimating basal area and stems per hectare using aerial auxiliary variables. Forest Science, 51(2), 109-119.
Liang, Xinlian & Hyyppä, Juha & Matikainen, L. (2007). Deciduous-coniferous tree classification using difference between first and last pulse laser signatures. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 36(3), 253-257.
McRoberts, R. E., Nelson, M. D. & Wendt, D. G. (2002). Stratified estimation of forest area using satellite imagery, inventory data, and the k-nearest neighbours technique. Remote Sensing of Environment, 82, 457-468.
Maltamo, M., Eerikäinen, K., Pitkänen, J., Hyyppä, J., & Vehmas, M. (2004). Estimation of timber volume and stem density based on scanning laser altimetry and expected tree size distribution functions. Remote Sensing of Environment, 90(3), 319-330.
10.1016/j.rse.2004.01.006
Maltamo, M., Eerikäinen, K., Packalén, P., & Hyyppä, J. (2006). Estimation of stem volume using laser scanning based canopy height metrics. Forestry, 79(2), 217-229.
10.1093/forestry/cpl007
