A Comparison of Linear-Mode and Single-Photon Airborne LiDAR in Species-Specific Forest Inventories

Single-photon airborne light detection and ranging (LiDAR) systems provide high-density data from high flight altitudes. We compared single-photon and linear-mode airborne LiDAR for the prediction of species-specific volumes in boreal coniferous-dominated forests. The LiDAR data sets were acquired at different flight altitudes using Leica SPL100 (single-photon, 17 points <inline-formula> <tex-math notation="LaTeX">$\cdot ~\text{m}^{-2}$ </tex-math></inline-formula>), Riegl VQ-1560i (linear-mode, 11 points <inline-formula> <tex-math notation="LaTeX">$\cdot ~\text{m}^{-2}$ </tex-math></inline-formula>), and Leica ALS60 (linear-mode, 0.6 points <inline-formula> <tex-math notation="LaTeX">$\cdot ~\text{m}^{-2}$ </tex-math></inline-formula>) LiDAR systems. Volumes were predicted at the plot-level using Gaussian process regression with predictor variables extracted from the LiDAR data sets and aerial images. Our findings showed that the Leica SPL100 produced a greater mean root-mean-squared error (RMSE) value (41.7 <inline-formula> <tex-math notation="LaTeX">$\text{m}^{3}~\cdot $ </tex-math></inline-formula> ha<sup>−1</sup>) than the Leica ALS60 (39.3 <inline-formula> <tex-math notation="LaTeX">$\text{m}^{3}~\cdot $ </tex-math></inline-formula> ha<sup>−1</sup>) in the prediction of species-specific volumes. Correspondingly, the Riegl VQ-1560i (mean RMSE = 33.0 <inline-formula> <tex-math notation="LaTeX">$\text{m}^{3}~\cdot $ </tex-math></inline-formula> ha<sup>−1</sup>) outperformed both the Leica ALS60 and the Leica SPL100. We found that the cumulative distributions of the first echo heights >1.3 m were rather similar among the data sets, whereas the last echo distributions showed larger differences. We conclude that the Leica SPL100 data set is suitable for area-based LiDAR inventory by tree species although the prediction errors are greater than with data obtained using the modern linear-mode LiDAR, such as Riegl VQ-1560i.

and backscattered signals. Discretized data are generated from waveforms, and the discretization of waveforms is usually implemented as waveform processing during data acquisition [4]. Discretized LiDAR data are easier to process and require less storage capacity than waveform data, which may be the reason why full-waveform data have been rarely applied in operational forest inventories. Typical linear-mode airborne LiDAR systems use a single wavelength and can record several echoes per emitted pulse. Ranging is based on the time elapsed between the emission of the pulse and the capture of the returning photon surge. Typical linear-mode LiDAR systems are also capable of recording intensity values, which are related to the strength of the backscattered signal. The amplitude of the backscattered signal is usually recorded as an intensity value [5]. The intensity value of the returning echo depends on the emitted LiDAR power, the range between the sensor and target, and the shape and directional reflectance of the target [6]. The intensity values can be used to discriminate certain tree species [7]- [10].
Single-photon (or photon counting) LiDAR systems are more sensitive in detecting backscattered photons than conventional linear-mode LiDAR systems [11]. Typically, a few photons are sufficient to record a range measurement. The greater sensitivity of single-photon LiDAR enables the collection of high-density point clouds from altitudes higher than that obtained with linear-mode LiDAR. For forest inventory purposes, linear-mode airborne LiDAR systems generally operate at altitudes of 500-2000 m [12], while single-photon LiDAR systems can operate at 4500 m above ground level [13]. As of 2020, few commercial single-photon airborne LiDAR systems are available.
The Leica SPL100 single-photon system (used in this study) splits the laser beam into an array of 10 × 10 collimated subbeams (beamlets) using a diffractive optical element and measures the backscattered signal for each beamlet [14]. In addition to beam splitting, the Leica SPL100 system utilizes extremely short laser pulses and has the capability to emit six million pulses per second. As a reference, a linear-mode Riegl VQ-1560i can carry out 1.3 million measurements per second. In practice, single-photon LiDARs enable the acquisition of a significantly greater number of measurements per unit area than linear-mode systems when flown at a similar altitude using the same flight speed. Since the sensor can record echoes based on just a few backscattered photons and the technical realization of the receiver is different from that of linear-mode LiDARs, the intensity measurements are not This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ directly comparable to linear-mode LiDAR systems. The Leica SPL100 system measures the intensity values based on the width of the returning echoes [15], [16].
The SPL100 operates in the visible green region (532 nm) of the spectrum that is suboptimal for the remote sensing of vegetation. In the green wavelength, the reflectance from vegetation is considerably weaker compared to near-infrared wavelengths [17]. The linear-mode LiDAR systems commonly operate at near-infrared wavelengths (e.g., 1064 or 1550 nm). For example, Kukkonen et al. [10] observed that, in forested areas, the multispectral linear-mode Optech Titan sensor has a considerably lower capability to capture range measurements in the green region than at near-infrared wavelengths. The Leica SPL100 can detect echoes with lower energy than conventional LiDAR systems, and therefore, the low reflectivity in vegetation is not as problematic with single-photon instruments compared to conventional linear-mode LiDAR systems. The green wavelengths are also more sensitive to solar noise than near-infrared wavelengths [11], [17], [18]. The detection sensitivity and green wavelength of the SPL100 system mean that a significant fraction of solar and atmospheric noise photons are recorded. Several algorithms have been developed to deal with this noise [17], [19], [20]. In practice, noise removal is carried out using the software provided by the sensor manufacturer.
It is assumed that single-photon LiDAR systems provide cost savings compared to linear-mode LiDARs, especially in large area data acquisitions. Typically, the flying altitude of single-photon LiDAR campaigns is considerably higher than that of linear-mode LiDAR systems. The flying speed with single-photon LiDAR can also be increased, compared to linear-mode sensors, because of the increased sensitivity and point density. For many technical and economic reasons, single-photon airborne LiDARs are still rare, and therefore, linear-mode LiDAR is typically preferred to single-photon LiDAR in practical campaigns so far.
Single-photon LiDAR data have been used in few studies in a forestry context. Swatantran et al. [17] investigated the capability of single-photon LiDAR data to characterize canopy structures and compared the results with those of linear-mode LiDAR. They observed that single-photon and linear-mode LiDAR data are able to characterize the canopy surface similarly, but that the pulses of linear-mode LiDAR penetrate deeper inside the canopy. Wästlund et al. [21] and Yu et al. [22] studied the estimation of forest variables using the Leica SPL100 and found that selected forest attributes can be predicted with comparable error rates with both linear-mode LiDAR and single-photon LiDAR data. Wästlund et al. [21] also found that single-photon LiDAR beams penetrate better through the canopy in hemiboreal forests than linear-mode LiDAR beams. Wästlund et al. [21] used the Leica SPL100 instrument, whereas Swatantran et al. [17] applied an experimental airborne single-photon LiDAR known as the High-Resolution Quantum LiDAR System. Mandlburger et al. [23] reported that linear-mode LiDAR data can better characterize the vegetation below the dominant canopy surface than Leica SPL100 data. Brown et al. [24] also reported that Leica SPL100 data poorly characterize forest canopy structures compared to linear-mode LiDAR. The studies emphasize that investigations of noise filtering in forest conditions are still needed. The noise removal algorithms used with single-photon LiDAR data may also eliminate vegetation echoes below the top surface of a forest canopy [23].
To date, no studies have evaluated the potential of single-photon airborne LiDAR data for the prediction of tree species-specific forest attributes although forest attributes are often needed by tree species. In general, aerial images are used in combination with airborne LiDAR data to incorporate more tree species information. Aerial image features are especially useful in separating coniferous and broadleaved tree species because they exhibit different reflectance at the near-infrared region of the spectrum [25]. Previous studies have also reported that intensity features extracted from LiDAR data deviate between tree species (e.g., [26] and [7]). It must be noted that single-photon LiDAR technology cannot measure intensity values based on the echo waveform amplitude, as is the case with linear-mode LiDAR. We refer to the intensity of the Leica SPL100 sensor as pseudointensity in this article to make a clear difference between linear-mode and single-photon LiDAR intensities. Features derived from aerial images usually outperform LiDAR intensity features in the prediction of forest attributes by tree species [27], [28]. When tree species are not considered, LiDAR-based features alone usually provide a sufficiently low prediction error for forest attributes.
In this study, our objective was to compare single-photon and conventional linear-mode airborne LiDAR in the prediction of species-specific forest attributes using an area-based approach (ABA). We used three LiDAR data sets: low-density linear-mode LiDAR data (Leica ALS60), high-density linearmode LiDAR data (Riegl VQ-1560i), and high-density singlephoton LiDAR data (Leica SPL100). The LiDAR data sets were compared in typical acquisition modes applied in forest inventories in Finland. This means that the Leica ALS60, which is an older-generation device, had substantially lower pulse density and was operated at a higher altitude than the newer generation Riegl VQ-1560i. The data acquisition parameters of Leica SPL100 were chosen, as the flight altitude of single-photon LiDARs is typically higher than linear-mode LiDARs, to exploit its capability to efficiently collect data from large areas. First, we investigated the differences in the distributions of echo heights and the differences in the intensity (pseudointensity in the case of SPL100) distributions in pure Scots pine (Pinus sylvestris [L.]), Norway spruce (Picea abies [L.] Karst.), and broadleaved plots. Second, we evaluated the performance of different LiDAR data sets in the prediction of plot volumes by tree species while coupling LiDAR features with the spectral features of aerial images.

A. Study Area and Sample Plot Data
The study area is located near the city of Tampere (N 61 • 10 E 23 • 52 ) in southern Finland (see Fig. 1). The study area is part of an operational forest inventory area where the Finnish Forest Centre had previously measured sample plots in the summer of 2017. The locations of the plots were determined using systematic stratified sampling [29]. The sample plots represent managed boreal forests with young, middle-aged, and mature development classes. The dominant tree species by timber volume in the plots are Norway spruce (47%) and Scots pine (37%). The remainder of the plots is dominated by broadleaved species, mostly silver birch (Betula pendula Roth.) or downy birch (Betula pubescens Ehrh.). Typically, broadleaved species also exist as mixtures in plots dominated by coniferous tree species.
The field data comprised 310 circular plots with either 9-or 12.62-m radius. We used the 12.62-m radius if the number of trees within the 9-m radius was <20. The 9-m radius was used in 70% of plots. The plots were accurately positioned using a GNSS system and an external antenna elevated to 5 m. We postcorrected the plot locations using a differential GNSS algorithm and reference positions. Diameter at breast height (DBH) and species were measured for each tree with a DBH ≥ 5 cm. The height of one sample tree (mean by DBH) of each species in each story class was measured in each plot. Heights for trees without height measurement were predicted using the species-specific height-DBH curve (mixed-effect model) proposed by Eerikäinen [30]. The model predictions were calibrated (i.e., localized) using the random parameters of the mixed-effect model [31]. We employed volume functions presented in [32] in order to calculate species-specific stem volumes for each tree using DBH and height as predictor variables. Volume (m 3 · ha −1 ) and basal area (m 2 · ha −1 ) per hectare were computed by tree species for each plot. The dominant tree species of each plot was determined according to species-specific volumes. A summary of the plot attributes is presented in Table I.

B. Remotely Sensed Data
We used two linear-mode LiDAR data sets and one single-photon LiDAR data set in this study. All data sets were collected under leaf-on conditions from a fixed-wing aircraft, either in 2017 or 2018. The first linear-mode LiDAR data set was collected with a Leica ALS60 scanner. According to the plot-level mean pulse densities (see Table II), it had a low pulse density (< 1 pts · m −2 ), which is commonly used in operational forest inventories in Finland [2]. The second linear-mode LiDAR data set was collected with a Riegl VQ-1560i scanner and had a greater pulse density (10.95 pts · m −2 ). The single-photon LiDAR data set was collected with a Leica SPL100 sensor and had the greatest pulse density (17.00 pts · m −2 ). It should be noted that the SPL100 sensor uses a nutating mirror (Palmer) scanning mechanism, which means that the point density is greater on the edges of the flight line than in the middle of the flight line. The Riegl VQ-1560i and Leica ALS60 apply rotating polygon and oscillating mirror scanning mechanisms, respectively. The scanning patterns of the LiDAR sensors and the LiDAR point clouds are illustrated in Fig. 2. The abbreviations hereafter used for the LiDAR systems and the detailed acquisition parameters are presented in the columns of Table II.
The data provider preprocessed the SPL100 data with Leica's own software tailored to this instrument, and this included noise removal. Leica Geosystems has not disclosed the technical details of their noise removal algorithm. We did not implement any additional noise removal steps.
The ALS60 and VQ1560i systems record multiple returns per pulse. The SPL100 system may also register multiple returns per emitted beamlet. The echo detectors of the SPL100 system operate similar to the detectors of linear-mode LiDARs, at least in a specific dynamic range [23]. The original echo categories were first-of-many, only, intermediate, and last-of-many. We added only echoes to both first-of-many and last-of-many echoes in order to establish first and last categories. The SPL100 system recorded only one echo per emitted beamlet more often than ALS60 or VQ1560i. The proportions of echoes by the original echo categories in the SPL100, VQ1560i, and ALS60 data sets are presented in Table III. The echoes were classified as ground and vegetation return separately for each LiDAR data set, following the approach proposed by [33]. The ground returns were used to  LIDAR DATA SETS, PLEASE REFER TO TABLE II interpolate a digital terrain model (DTM) using a Delaunay triangulation. The echoes backscattered from vegetation were height-normalized by subtracting the DTM from the echo heights.
The intensities of echoes are determined differently with the linear-mode and single-photon sensors, which must be taken into account in the comparison of intensity response. Linear-mode LiDAR systems (here ALS60 and VQ1560i) digitize discrete echoes from each returned waveform. Usually, the peak amplitudes of the returning waveform are recorded as discrete-return intensities, which indicates the strength of each discrete echo [5]. The intensities of ALS60 data are directly related to the peak amplitude of the received waveform, and an automatic gain control was not used, i.e., the gain was fixed in the acquisition of the ALS60 data. The algorithms applied in the digitization step are not available to end-users. We normalized the intensities of ALS60 data based on the range between the target and LiDAR instrument. The range-normalization was similar to Korpela et al. [6] and Gatziolis [34]. The intensity values were corrected using an exponent 2.0 and 2450-m mean elevation above ground. The SPL100 device cannot measure the peak amplitude of an echo due to the operating principle of the single-photon LiDAR system. However, the Leica SPL100 generates a pseudointensity using the estimated width of the backscattered pulse. The width of the pulse is estimated by measuring the elapsed time at which the signal (backscattered pulse) exceeds and drops below the detection threshold [15]. The intensities of the VQ1560i data were determined using an amplitude relative to the amplitude seen at the detection threshold. The intensity measurements of the VQ1560i instrument were processed using the "relative reflectance" algorithm described by Riegl [35].
We also used aerial images that are typically applied in species-specific forest inventories [2]. The aerial images used in this study were collected between June 30 and July 1, 2017. The images were captured at an altitude of 4000 m above ground level (longitudinal overlap 80%) using a Z/I Imaging (01-0128) digital aerial camera that is capable of recording red, green, blue, and near-infrared bands. We also computed the normalized difference vegetation index (NDVI) that is derived from the recordings at the red and near-infrared wavelengths [36]. The camera had a focal length of 30 mm and a 1920 × 3456 pixel sensor for color images. The ground sampling distance of the aerial images was 1.6 m for the color bands. The aerial images were not pansharpened or orthorectified because photogrammetric techniques were applied in the feature extraction [37].

C. LiDAR and Aerial Image Features
We extracted features from the LiDAR data sets for the sample plots using our own software implemented in R environment [38]. The features were computed by echo categories of first, last, and intermediate. The height-related features computed from the first and last echoes consisted of the following statistics: maximum, mean, median, standard deviation, percentiles (10%, 20%, . . ., 90%), skewness, and kurtosis. Moreover, canopy density features were computed at fixed heights of 0.5, 2, 5, 10, and 15 m. The aforementioned features were also computed with intensity values (first and last echoes), except for the minimum of intensity values and density features. The intermediate echoes were used to compute mean, standard deviation, percentiles (20%, 40%, . . ., 80%), and echo proportions. Hereafter, the LiDAR features, other than intensity features, are referred to as LiDAR structure features. A 1.3-m height cutoff was used in the computation of all features other than densities. Intermediate echoes were rare in the SPL100 data (only 0.1% of the echoes), and, therefore, we did not compute intermediate features for the SPL100 data.
Each LiDAR echo classified as the category first with a height above the plot-level median height was backprojected into the unrectified aerial images using external and internal orientations, and digital number (DN) values were retrieved from the images relative to the echoes [37]. The cutoff value of the plot-level median was used to remove the effects of ground and vegetation below the dominant tree layer. Each echo was linked to several images due to image overlap. Therefore, a mean DN value by the band was computed for each echo. Finally, plot-level means and standard deviations were computed from the mean DN values for the red, green, blue, and near-infrared bands and NDVI. The features extracted from the aerial images are hereafter referred to as optical image features.
All the extracted features are listed in Table IV. For the sake of clarity, we categorized the features into three groups: LiDAR intensity/pseudointensity, LiDAR structure, and optical image features. The LiDAR structure group comprises all the LiDAR features other than intensity features. All the features were computed separately and similarly for each LiDAR data set.

D. Echo Height and Intensity Distributions in Pure Species Plots
We examined the echo height and intensity/pseudointensity distributions in pure tree species plots, i.e., separately in spruce, pine, and broadleaved forests. This was done for each LiDAR data set. Height distributions were examined  IV   FEATURES EXTRACTED FROM THE LIDAR DATA SETS. THE FEATURES ARE CATEGORIZED INTO THREE GROUPS: LIDAR  INTENSITY/PSEUDOINTENSITY, LIDAR STRUCTURE, AND OPTICAL IMAGE FEATURES. THE COMBINED SET OF FEATURES  TAKEN FROM LIDAR INTENSITY/PSEUDOINTENSITY AND LIDAR STRUCTURE CATEGORIES IS REFERRED TO AS  "ALL LIDAR FEATURES." THE PSEUDOINTENSITY REFERS TO THE INTENSITY MEASUREMENTS OF THE LEICA  SPL100 INSTRUMENT. THE FEATURES WERE COMPUTED BY ECHO CATEGORIES (SEE SECTION II-C). ABBREVIATIONS:  I-INTENSITY, H-HEIGHT, C-ECHO COUNT, OI-OPTICAL IMAGE, R-RED, G-GREEN, B-BLUE,  NIR-NEAR-INFRARED, AND NDVI-NORMALIZED DIFFERENCE VEGETATION INDEX using cumulative distributions, which indicates how the pulses penetrate through the canopy layers and how the heights of LiDAR echoes accumulate. We left intermediate echoes out of the analysis because they were absent in many plots in the case of the SPL100 and ALS60 data sets. The intensity/pseudointensity distributions were compared to examine the differences by tree species. This indicates the potential of intensity/pseudointensity to discriminate between tree species.

E. Estimation Methodology
We predicted the species-specific volumes using the multivariate Gaussian process (GP) regression [39]. The GP regression is a widely used machine learning method that is related to kriging and support vector machines. In the GP regression, the relationship between the predictor variables and the response variables is modeled as a stochastic process described by a covariance function. The main benefits of the GP regression are its relative ease of use, computational tractability, and attractive theoretical properties; for a detailed description, see [40]. The GP regression has provided promising results in similar settings to those in this study, with prediction errors even lower than the nearest neighbor approach tailored for the Finnish species-specific forest inventories [2]. In this GP variant, several response variables could be modeled simultaneously. The predicted total volume was computed as a sum of the predicted species-specific volumes. The postprocessing of negative predictions was implemented, as presented in detail by Varvia et al. [39].
The covariance function is a key component of the GP regression approach since it ultimately defines the function to be trained and its properties, such as smoothness and underfitting or overfitting. We applied a stationary Matern covariance function that can be expressed with ν = 2/3 and σ = 1, as follows: where d is the Euclidean distance between data points and l is a correlation length hyperparameter. The GP regression also has an error variance parameter e that determines the magnitude of noise in the training data. We optimized the hyperparameters e and l associated with the GP regression using a grid search algorithm, which was executed following the leave-plot-out principle. The mean of root-mean-squared errors [RMSE; see (2)] value associated with the predicted species-specific volumes was a loss function. We carried out the grid search separately for each LiDAR data set and set of predictor variables. The candidate values for e and l in the grid search were selected from a predefined feasible range to avoid overfitting. The candidate values for e were from the range [0.08, 0.18] using an interval of 0.02. Correspondingly, the candidate values for l were from the range [4,18] using an interval of 1. Altogether, the grid search tested 84 combinations of hyperparameters for each data set.
We performed the prediction of species-specific volumes using the following sets of predictor variables: 1) LiDAR structure; 2) all LiDAR features (structure + intensity/pseudointensity features); and 3) all LiDAR and optical image features. Due to the moderate number of predictor variables, we did not use variable selection routines, i.e., all the features of a set were always used. Varvia et al. [39] used the GP regression in the prediction of species-specific forest attributes and did not use variable selection either.

F. Accuracy Assessment
We evaluated the predictive performance of the GP regression using the predictions obtained through leave-one-plotout cross validation. The predictive performance of the GP regression was assessed using RMSE and mean difference (MD) whereŷ i is the predicted value in plot i , y i is the observed value in plot i , and n is the number of plots.

A. Echo Height and Intensity Distributions
We found that the cumulative distributions of echo heights (above 1.3 m) were dissimilar between the LiDAR data sets. In pure spruce and pine forests, the cumulative distributions of linear-mode first echo heights are almost similar compared to that of SPL100 (see Fig. 3). In pure broadleaved forests, the SPL100 sensor detected a larger proportion of first echoes from low vegetation (<5 m) than the linear-mode sensors. In the case of the last echoes, the distribution curves associated with the ALS60 data set were dissimilar to those of the SPL100 and VQ1560i data sets. The ALS60 sensor detected a larger proportion of last echoes from low vegetation than the SPL100 and VQ1560i sensors regardless of tree species.
We also report the proportion of first and last echoes ≤ 1.3 m in Table V. The findings showed that the SPL100 data set had a larger proportion of echoes returning from heights ≤ 1.3 m in pine and broadleaved forests than ALS60 and VQ1560i.
A closer look at the intensity (pseudointensity in the case of SPL100 data) distributions revealed that all three LiDAR systems provided unequal intensity distributions in the pure species plots (see Fig. 4). In particular, the linear-mode VQ1560i and ALS60 behaved differently, as the VQ1560i provided greater intensities for the spruce plots compared to pine, while the situation was generally the opposite for ALS60. SPL100 provided similar pseudointensity distributions for coniferous species. However, the pseudointensity distribution of SPL100 in the broadleaved plots was different compared to the coniferous plots. The mean of pseudointensities was generally smaller in the broadleaved plots than in the coniferous plots.

B. Volume Prediction by Tree Species
Figs. 5 and 6 show the RMSE and MD values associated with the species-specific and total volume predictions using the following feature sets: LiDAR structure, all LiDAR features, and all LiDAR features and optical image features. In Section III-A, we showed that some intensity distributions were dissimilar between tree species, and the differences in the intensity distributions were most visible in the VQ1560 data in the coniferous plots. The results of the volume predictions also indicated that intensity features are beneficial in the prediction of volume by tree species, especially with the linear-mode LiDARs. The benefit of the inclusion of intensity (pseudointensity in the case of SPL100 data) features was greater with the ALS60 and VQ1560i data than with the SPL100 data. The findings indicate that the pseudointensity of SPL100 was not useful in the prediction of species-specific volumes. This finding was also in line with our preliminary results regarding the classification of dominant tree species using the LiDAR intensity features (results not shown here). With LiDAR structure features only, the means of species-specific RMSE by LiDAR data sets were 51.5, 46.0, and 57.8 m 3 · ha −1 for SPL100, VQ1560i, and ALS60, respectively. Corresponding figures with all LiDAR features were 50.3, 40.3, and 51.6 m 3 · ha −1 for SPL100, VQ1560i, and ALS60.
The results also indicate that LiDAR features, other than intensity, may have a role in the separation of tree species when volumes are predicted by tree species (see high-density VQ1560i and low-density ALS60 data on the upper row in Fig. 5). The MD values associated with the species-specific predictions were not particularly high in any case. The MD for total volume increased when the number of predictor variables increased (see Fig. 6).
The combined use of all LiDAR and optical image features provided substantially smaller error rates for species-specific predictions than the use of LiDAR features only. The LiDAR sensor also had a clear effect on predictions in this case. The means of species-specific RMSE values by LiDAR data sets were 41.7, 33.0, and 39.3 m 3 · ha −1 for SPL100, VQ1560i, and ALS60, respectively, when using all LiDAR and optical image features. VQ1560i was superior to the other two instruments, in particular with respect to the RMSE values of pine and spruce, while SPL100 produced the largest RMSE value with respect to all attributes. With respect to total volume, the combined use of all LiDAR and optical image features performed similarly for all of the instruments. The RMSE value associated with total volume predictions decreased slightly in every instance when optical image features were added to the set of predictors. The MD value slightly increased when optical image features were added to the set of predictors. Predicted versus observed volumes by tree species and LiDAR instruments (when all LiDAR and optical image features were used) are shown in Fig. 7.

IV. DISCUSSION
The VQ1560i LiDAR data outperformed the ALS60 and SPL100 data in the prediction of species-specific volumes with the data acquisition settings used in this study. The use of  Table II. OI-optical image. linear-mode LiDAR intensities notably decreased the prediction errors of species-specific plot volumes compared to the use of LiDAR structure features only. The pseudointensities of the SPL100 data set are not entirely comparable to the intensities provided by linear-mode LiDARs, but we see the comparison as relevant because intensities and pseudointensities may be used in a similar manner. The results showed that the pseudointensities are not very useful in the prediction of species-specific plot volumes. However, we found that the pseudointensity distributions of broadleaved and coniferous plots were distinct, which means that the pseudointensities might be beneficial in the discrimination of coniferous and broadleaved tree species. The benefit of intensity features may, however, be limited in operational inventories where optical image features are also used [2]. As our results showed, the combined use of LiDAR and optical image features clearly provides the lowest RMSE values. Linear-mode VQ1560i data provided the smallest RMSE values in the prediction of species-specific volumes when both LiDAR and optical image features were used. Single-photon SPL100 data provided the greatest RMSE values, and the linear-mode ALS60 data slightly outperformed it. In general, the prediction errors of species-specific volumes were consistent with previous studies conducted in similar forests [39], [41].
The LiDAR structure features had an effect on the predictive performance as well, which indicates that the data sets characterized the forest canopies somewhat differently. With the data acquisition settings used in this study, VQ1560i clearly outperformed ALS60 and slightly outperformed SPL100 when species-specific volumes were predicted using LiDAR structure features only. The distributions of echo heights collected from forests using single-photon and linear-mode LiDAR sensors have been compared previously [22]- [24]. In the previous studies, it has been observed that the pulses of linear-mode LiDAR penetrate more inside the forest canopy than the pulses of single-photon LiDAR. Yu et al. [22] showed that tree species composition, the vertical structure of the forest, and the flight altitude of data acquisition have a considerable impact on the distribution of echo heights when a single-photon LiDAR system is used. The results of our study suggest that the differences in the cumulative distributions of the first echo heights (considering echo heights above 1.3 m) between linear-mode and single-photon sensors are minor in boreal forests. In the case of the last echoes, the cumulative distribution curves of echo heights associated with ALS60 were steeper at low heights than with the other sensors, especially in pine and broadleaved forests (see Fig. 3). This could be possibly due to the large footprint of ALS60. ALS60 did not record as large a proportion of last echoes from the upper canopy than the other instruments, and therefore, echoes from low heights are emphasized in the cumulative distribution. Note that we did not analyze cumulative distribu- Fig. 6. MD (m 3 · ha −1 ) values associated with the volumes predicted by tree species and total volume. Row indicates the features used, and column indicates the LiDAR system used. For the abbreviations of LiDAR systems, please refer to Table II. OI-optical image.
tions of intermediate echo heights since they were rare in the SPL100 and ALS60 data sets and were not particularly useful in the prediction of volumes (results not shown).
We found that the proportion of echo heights ≤ 1.3 m was greater in the SPL100 data set compared to the ALS60 and VQ1560i data sets in pine and broadleaved forests. A closer investigation of the SPL100 data reveals that the number of echoes backscattered from heights of 0-1.3 m above the ground level was large. Most of those near-ground echoes, correctly, were not classified as ground echoes. A large proportion of near-ground echoes may be due to noise in the characterization of the planar surface, which can be explained by the high sensitivity of single-photon sensors [23], [24]. Wästlund et al. [21], Yu et al. [22], and Brown et al. [24] also reported that single-photon LiDAR data sets have more echoes from near ground level than linear mode LiDAR data. Brown et al. [24] suggested that a large number of emitted echoes in the SPL100 data is the main reason for the large number of echoes backscattered from the ground (echoes per unit area). Mandlburger et al. [23] compared single-photon and linear-mode LiDAR data sets with similar point densities. The results reported by Brown et al. [24] concerning the number of ground echoes were in agreement with Mandlburger et al. [23]. In our SPL100 data set, one reason for the large proportion of near-ground echoes might be noise removal. Our hypothesis is that noise removal designated many echoes from the canopy and between the ground and canopy layer as noise echoes although the echoes near the ground were not identified as noise echoes. This increases the proportion of near-ground echoes, which should be taken into account in the calculation of features for the modeling of forest attributes. In this study, we used the height cutoff at 1.3 m in the extraction of features.
LiDAR intensities are used in species-specific forest inventories although they are generally outperformed by aerial images in terms of predictive power in the estimation of species-specific forest attributes. We suggest that the differences in the LiDAR intensity responses are a major reason for the success of VQ1560i in the prediction of species-specific volumes using all LIDAR features. LiDAR intensities are commonly studied in the classification of tree species [7], [10]. We found that the intensity distributions of the VQ1560i data were distinct, particularly in spruce and pine forests. The shapes of the intensity distributions were also very different in spruce and pine forests compared to the other linear-mode instrument (ALS60). We suggest that the reasons for the differences between spruce and pine forests with linear-mode LiDARs are related to intensity processing, data acquisition parameters, and tree phenology. In general, the linear-mode instruments VQ1560i and ALS60 provided differently shaped intensity distributions. In the SPL100 data, the differences in the pseudointensity distributions in the coniferous and broadleaved forests can be explained by how the pseudointensity is derived from the echo width. It would appear that the approximated echo width is shorter in the broadleaved forest than in coniferous forests. Broadleaved forests typically have relatively solid and dense canopy surfaces covered by leaves, whereas the needles of coniferous species are organized hierarchically into shoots that display more volumetric reflectance characteristics [42].
Sensor manufacturers use different algorithms to determine LiDAR intensity from waveform data [43]. Therefore, even among linear-mode LiDARs, the intensities of different sensor brands are rarely fully comparable. The intensity is often determined for each detected echo by the amplitude of the echo signal although the amplitude can be determined in different ways. Additional processing steps, such as amplitude normalizations or corrections, may be implemented in the data processing routines. These processing steps are rarely explained in detail to end-users. Riegl [35] has published a manual detailing the processing routines that are also used for VQ1560i data. The calibrated amplitude and relative reflectance associated with the echo signals are introduced in the manual. In our study, the intensities of the VQ1560i data were determined following the principles of relative reflectance described in the manual. The determination of the relative reflectance may partly explain the success related to the separation of pine-and spruce-dominated plots reported here. Determination of the relative reflectance aims to solve the issues related to range dependence and directionality of the reflections [35]. Therefore, the relative reflectance readings facilitate the interpretation of intensity measurements collected from different scan positions. We used range-normalization for the intensities of the ALS60 data set in order to make the comparison between linear-mode instruments as fair as possible.
In addition to intensity processing, the superiority of VQ1560i compared to ALS60 in the separation of tree species may be related to the sensor and data acquisition parameters. The VQ1560i data were collected from a lower altitude than the ALS60 data, resulting in a smaller footprint (36 cm versus 56 cm, respectively) and considerably greater density (10.95 versus 0.63 pts · m −2 , respectively). Previous studies have indicated that area-based predictions are rather independent to pulse density [10], [44], [45]. However, a lower flying altitude produces a narrower footprint of the laser pulse at the target, which affects the features extracted from the LiDAR data [46]. In general, the footprint size affects how deeply the laser beams are able to penetrate the tree crowns [47]. Penetration also varies according to tree species; for example, spruce canopies are more difficult to penetrate than pine canopies. The smaller footprint associated with the VQ1560i data (see Table II) may partly explain the differences in the intensity distributions in the pure spruce and pine plots (see Fig. 4).
In this study, we did not compare different LiDAR sensors using similar data acquisition parameters, e.g., using the same flying altitude. This must be kept in mind when the results are interpreted. We aimed to compare LiDAR sensors using the acquisition parameters that have been or would be used with the sensor in question in the operational remote sensing-based forest inventory in Finland. For example, if the ALS60 data would have been acquired from the altitude that produced a similar point density as the VQ1560i data, the intensity distributions and penetration through the upper canopy layer could probably have been very different.
The ALS60 data were acquired at the end of June 2017, and the VQ1560i and SPL100 data sets were acquired during a narrow time period at the end of May 2018. All data sets were collected under full leaf-on conditions. However, the phenology of coniferous tree species may cause different intensity responses in early and mid-summer. Hovi et al. [48] explored tree-level intensity values extracted from full-waveform airborne LiDAR data. Their results indicated that the distributions of intensity measurement between pine and spruce are more distinct in early than late summer. These phenological effects may partly explain the different intensity distributions between spruce and pine with the VQ1560i and ALS60 systems (see Fig. 4).
Single-photon LiDAR sensors are targeted at applications that require a high point density over large areas. Previous studies, however, indicate that increased point density does not usually decrease prediction errors in forest inventories that rely on ABA, at least when tree species are ignored [10], [22], [44], [45]. The resulting point density depends on acquisition parameters, such as flying altitude, which may affect the characteristics of LiDAR measurements. The values of the data acquisition parameters can significantly affect the prediction errors of the forest attributes, as has been noticed with linear-mode instruments [12], [49]. This study applied ABA in the prediction of species-specific volumes. The applicability of single-photon LiDAR data in tree-level species-specific inventories should also be evaluated.

V. CONCLUSION
Both linear-mode and single-photon LiDAR data sets characterize forest canopies sufficiently for the purposes of ABA inventories. The intensity distributions of Riegl VQ-1560i data were different from those of Leica ALS60 data, which would indicate that Riegl VQ-1560i is able to distinguish between pine and spruce forests. The pseudointensity distributions of Leica SPL100 were dissimilar in broadleaved and coniferous forests. This study examined the operational applicability of the LiDAR sensors in Finland and showed that the linear-mode Riegl VQ-1560i (high-density) outperformed the Leica SPL100 (high-density) and Leica ALS60 (lowdensity) in the prediction of tree species-specific volumes, both with and without the optical image features. The Riegl VQ-1560i sensor provided smaller prediction errors for pine and spruce volumes than the other sensors. In general, the Leica SPL100 provided the largest prediction error associated with species-specific volumes. Noise removal of single-photon data is critical and may remove the LiDAR pulses backscattered from the canopy, which are not true noise echoes. Our findings indicate that single-photon LiDAR data could substitute the linear-mode LiDAR data in operational species-specific forest inventories, provided that the minor increases in prediction errors are acceptable. The findings showed that the combination of linear-mode LiDAR data and aerial images provides the smallest error rates associated with the species-specific timber volumes. Petri Varvia received the M.Sc. and Ph.D. degrees in applied physics from the University of Eastern Finland, Joensuu, Finland, in 2013 and 2018, respectively.
He was a Post-Doctoral Researcher with the Unit of Computational Sciences, Tampere University, Tampere, Finland. He is a Post-Doctoral Researcher with the School of Forest Sciences, University of Eastern Finland. His scholarly interests include statistical inverse problems, Bayesian statistics, and remote sensing.
Lauri Korhonen received the D.Sc. degree in forestry from the University of Eastern Finland (UEF), Joensuu, Finland, in 2011.
He is a Senior Researcher of forest mensuration with the School of Forest Sciences, UEF. His research interests include the measurements of forest canopy properties and forest inventory with airborne and spaceborne light detection and ranging (LiDAR) data.
Pekka Savolainen studied mathematics and physics at the University of Turku, Turku, Finland, and the Helsinki University of Technology (today Aalto University), Espoo, Finland.
He has been working with remote sensing and forest inventory for 30 years. He is an Research and Development Manager of Terratec Oy, Helsinki, taking part in forest inventories and data acquisition. In recent years, he has been occupied in developing a new reference measurement method, which registers precise positions for all tallied trees.
Matti Maltamo was born in Jyväskylä, Finland, in 1965. He received the M.Sc., Lic.Sc., and D.Sc. degrees (Hons.) in forestry from the University of Joensuu, Joensuu, Finland, in 1988Finland, in , 1992Finland, in , and 1998 He is a Professor of forest mensuration science with the Faculty of Science and Forestry, University of Eastern Finland, Joensuu. He has also been a Visiting Professor with the research group of Prof. Erik Naesset at the Norwegian University of Life Sciences, Ås, Norway. He was together with Naesset and prof Jari Vauhkonen, the editor of the textbook, Forestry Applications of Airborne Laser Scanning-Concepts and Case Studies (Springer, 2014). He has authored over 200 scientifically refereed articles. His research interests include forestry applications of airborne laser scanning (ALS).
Dr He is the Professor of optimization of forest management with the Faculty of Science and Forestry, School of Forest Sciences, University of Eastern Finland, Joensuu. Since 2007, he has been a consultant for remote sensing-based forest inventory. He has authored more than 100 peer-reviewed research articles. His research interests include both practical and theoretical aspects of utilizing remote-sensing data in the monitoring and assessment of the forest environment.