Comparison of Grid-Based and Segment-Based Estimation of Forest Attributes Using Airborne Laser Scanning and Digital Aerial Imagery

Abstract

Forest management planning in Finland is currently adopting a new-generation forest inventory method, which is based on interpretation of airborne laser scanning data and digital aerial images. The inventory method is based on a systematic grid, where the grid elements serve as inventory units, for which the laser and aerial image data are extracted and the forest variables estimated. As an alternative or a complement to the grid elements, image segments can be used as inventory units. The image segments are particularly useful as the basis for generation of the silvicultural treatment and cutting units since their boundaries should follow the actual stand borders, whereas when using grid elements it is typical that some of them cover parts of several forest stands. The proportion of the so-called mixed cells depends on the size of the grid elements and the average size and shape of the stands. In this study, we carried out automatic segmentation of two study areas on the basis of laser and aerial image data with a view to delineating micro-stands that are homogeneous in relation to their forest attributes. Further, we extracted laser and aerial image features for both systematic grid elements and segments. For both units, the feature set used for estimating the forest attributes was selected by means of a genetic algorithm. Of the features selected, the majority (61–79%) were based on the airborne laser scanning data. Despite the theoretical advantages of the image segments, the laser and aerial features extracted from grid elements seem to work better than features extracted from image segments in estimation of forest attributes. We conclude that estimation should be carried out at grid level with an area-specific combination of features and estimates for image segments to be derived on the basis of the grid-level estimates.

1. Introduction

In Finland, the forest inventory for forest management planning has traditionally been based on visual inventory by stands. In this method, the forest stands that are delineated on the basis of aerial photographs and their growing stock and site-related characteristics are measured or estimated in the field. The method requires a large amount of fieldwork. Therefore, the visual inventory method is to be replaced with a new-generation forest inventory method.

The new-generation forest inventory method will be based on interpretation of airborne laser scanning (ALS) data and digital aerial imagery using field sample plots as reference data. Laser scanning has been considered the most promising remote sensing technology in forest inventory, and it has been widely applied for stand-level forest inventories (e.g., [1,2,3,4]). On the other hand, ALS data are not well suited to estimation of tree species proportions or dominance at the applied point density (e.g., [5]). Accordingly, optical imagery is needed to complement the ALS data. Spectral features of optical images are typically used for separating different tree species, whereas textural features of optical images are mainly connected with the size and spatial arrangement of the tree crowns. In Finland, aerial images have been widely used in forest inventory since the 1950s [6], and their affordability and availability are good (e.g., [3,7]).

Statistically, the new-generation forest inventory method is based on two-phase sampling with stratification, where the inventory database is based on a systematic grid of sample units (i.e., grid elements as sample units), and the size of the grid elements should correspond to the size of field plots. Field measurements are allocated into strata that are commonly derived on the basis of earlier stand inventory data. Typical remote sensing data sources used in the new-generation forest inventory system are low-density ALS data (typically 1–2 pulses/m²) and digital aerial imagery with a spatial resolution of approximately 0.5 m containing the following spectral channels: blue (B), green (G), red (R) and near-infrared (NIR).

The forest variables are estimated for each element of the inventory grid, and each grid element is commonly defined as a square area. As an alternative to the grid-based approach, use of automatic stand delineation by image segmentation has been studied for defining inventory units (e.g., [8]). Automatic segmentation of ALS data has also been applied for delineating very small segments with a view to detecting individual tree crowns (e.g., [9,10]). In estimating stand level forest attributes (area based approach) the segment size is larger than single tree crowns. Stands delineated automatically on the basis of remote sensing imagery (i.e., image segments) have an advantage over grid elements (e.g., [11,12]). They can be delineated in such a way that they exactly follow the actual stand borders, whereas the grid elements are spatially ‘sparse’ in relation to the actual stand borders in the forest and they do not follow the borderlines accurately, instead often intersecting trees from more than one stand (e.g., [11]). On the other hand, the grid elements are unambiguously defined by their coordinates, so the same units can be used in subsequent inventories.

In delineation of the forest stands, the primary input variables are the mean height of the trees and the tree species composition (or dominance). From the stand delineation perspective, stand density usually is a secondary parameter. The height of the trees can be derived on the basis of the ALS data. As stated above, ALS data with the applied pulse density do not serve well the purpose of recognition of tree species. Therefore, again, optical aerial imagery is generally used for distinguishing forest stands on the basis of their tree species composition.

The geometrically three-dimensional nature of ALS data makes it possible to extract a large number of statistical features. When combining the ALS data with aerial photograph data, one finds that the number of available features increases further. Consequently, the dimensionality of the feature space increases greatly, and the data become sparse in relation to the feature space dimensions and the contrast between objects in the feature space weakens, making, for example, the nearest neighbor search unstable [13,14]. For the estimation procedure, the dimensionality of data must therefore be reduced, and a subset of features with good discrimination ability found. In an ideal case, the analyst would be able to infer the optimal feature combinations from the characteristics of the independent and the dependent data, but in the case of a physically complex and varying object, such as a forest, this is not possible, and automated feature selection methods must be used.

The objectives of this study were: (a) to find a suitable combination of laser and aerial data features for automatic stand delineation; (b) to find a suitable combination of laser and aerial data features for the estimation of forest attributes; and (c) to compare grid elements and automatically delineated stand polygons in the estimation of forest attributes.

2. Materials and Methods

2.1. Study Areas

The laser-scanning and aerial-image-based estimation was tested in two study areas. Study area 1 was located in the municipality of Lammi, in Southern Finland (approximately 61°19'N and 25°11'E). The area covered approximately 1,800 ha of state-owned forest. The field data in Study area 1 consisted of 281 fixed-radius (9.77 m) circular field sample plots that were measured in 2007. The plots were located with Trimble's GEOXM 2005 Global Positioning System (GPS) device, and the locations were processed with local base station data, resulting in an average error of approximately 0.6 m. Study area 2 was in Eastern Finland, in the municipalities of Kuopio and Karttula (approximately 62°55'N and 27°12'E), covering approximately 36,700 ha of mainly privately owned forest. The field data consisted of 546 fixed-radius (9 m) sample plots measured in 2009. In order to cover all types of forest, both study areas were stratified on the basis of earlier stand inventory data and the field sample plots were assigned to these strata. The location of the study areas and the sample plot layouts are presented in Figure 1.

Figure 1. Location of the study areas in Finland and of the sample plots within the study areas.

Figure 1. Location of the study areas in Finland and of the sample plots within the study areas.

Some differences were evident between the forest characteristics of the two study areas. In Study area 1, the total growing stock was more evenly distributed among the following tree species groups: Scots pine, Norway spruce, and deciduous trees, whereas Study area 2 was clearly dominated by Norway spruce. Furthermore, Study area 2 had a somewhat higher average stand volume, as well as greater variation in sample plot volumes. The statistics of the two study areas based on the sample plot measurements are presented in Table 1.

Table 1. Forest statistics of the study areas: average, maximum (Max.) and standard deviation (Std.) of the sample plot values.

	Study area 1			Study area 2
	Average	Max.	Std.	Average	Max.	Std.
Total volume, m3/ha	178.7	575.4	115.4	191.3	798.5	131.5
Volume of Scots pine, m3/ha	69.8	560.6	86.9	47.7	561.8	78.9
Volume of Norway spruce, m3/ha	63.7	575.4	94.9	102.9	739.2	128.0
Volume of deciduous species, m3/ha	45.2	312.0	56.2	40.7	400.4	63.9
Basal area, m2/ha	19.8	45.5	10.3	22.3	62.0	11.2
Mean height, m	17.0	30.5	6.7	16.9	35.6	6.7
Mean diameter, cm	21.1	50.2	9.4	20.7	60.3	10.0

Table 1. Forest statistics of the study areas: average, maximum (Max.) and standard deviation (Std.) of the sample plot values.

	Study area 1			Study area 2
	Average	Max.	Std.	Average	Max.	Std.
Total volume, m3/ha	178.7	575.4	115.4	191.3	798.5	131.5
Volume of Scots pine, m3/ha	69.8	560.6	86.9	47.7	561.8	78.9
Volume of Norway spruce, m3/ha	63.7	575.4	94.9	102.9	739.2	128.0
Volume of deciduous species, m3/ha	45.2	312.0	56.2	40.7	400.4	63.9
Basal area, m2/ha	19.8	45.5	10.3	22.3	62.0	11.2
Mean height, m	17.0	30.5	6.7	16.9	35.6	6.7
Mean diameter, cm	21.1	50.2	9.4	20.7	60.3	10.0

2.2. Remote Sensing Data

In Study area 1, the remote sensing data consisted of orthorectified color-infrared digital aerial imagery (containing near-infrared, red, and green bands) with a ground resolution of 0.5 m and ALS data acquired from a flying altitude of 1,900 m with a density of 1.8 returned pulses per square meter. In Study area 2, the remote sensing data consisted of orthorectified digital aerial imagery containing near-infrared, red, green, and blue bands with a ground resolution of 0.5 m and ALS data acquired from a flying altitude of 2,000 m with a density of 0.6 returned pulses per square meter.

In addition to use of the ALS point data, the ALS data were interpolated to a raster image format, for two output images: height and intensity. The raster-image pixel values of the height and intensity images were calculated with ArcGis Spatial Analyst tools, using inverse distance weighted (with a power of 2) interpolation based on the two nearest ALS points. The output laser images were resampled to a spatial resolution similar to that of the aerial images.

2.3. Automatic Image Segmentation

Stand delineation was carried out in the study areas via automatic segmentation of aerial images and ALS data interpolated to raster format. The segmentation was carried out in two phases. In the first phase, initial segmentation was performed via a modified implementation of the ‘segmentation with directed trees’ algorithm, which employs the local edge gradient [15,16]. The objective is to find all potential segment borders in this phase. Therefore, this method typically produces a very large number of small polygons when one is using high-resolution remote sensing data. In this study, the initial segmentation was based entirely on ALS height, corresponding mainly to stand height [17] (see Figure 2(a,b)). Prior to the initial segmentation, the ALS height data were pre-processed by Gaussian smoothing. The size of the smoothing window was 3 × 3 pixels, and five sequential smoothing operations were used, aimed at diminishing the within-stand variation and emphasizing between-stands variation.

In the second phase, the initial segments were processed via a region-merging algorithm that was guided by parameters such as the desired minimum size of the final segments and the similarity or dissimilarity of the segments to be merged [16]. The merging of regions into the final segments was carried out on the basis of laser height, laser intensity, and the NIR/R ratio of the aerial images, with the aim of taking into account also the tree species composition of the initial segments (Figure 2(c)). Two automatic segmentations, with minimum segment sizes of 350 m² and 0.1 ha, were carried out in both study areas.

Figure 2. Laser height data, initial segments based on the laser height data, and final segments (min. size 0.1 ha) based on a combination of laser and aerial image data. Examples from Study area 2.

Figure 2. Laser height data, initial segments based on the laser height data, and final segments (min. size 0.1 ha) based on a combination of laser and aerial image data. Examples from Study area 2.

2.4. Extraction of Laser and Aerial Image Features

Three remote sensing feature data sets were extracted for each of the study areas. In these sets, the remote sensing features were allocated to each sample plot from a square window or a segment in which the sample plot was located. The feature set Grid was extracted from a 20 × 20 m square window centered on each sample plot. The feature set Seg350 was extracted from image segments whose minimum size was set as 350 m². Feature set Seg1000 was extracted from image segments with a minimum size set as 0.1 ha.

The following statistical and textural features (max. 174) were extracted from the aerial images and ALS height and intensity of first pulse data for each feature data set:

• 1. Averages of pixel values of grid elements (20 × 20 m) and image segments surrounding each plot.
• 2. Standard deviations of pixel values of blocks, into which a 32 × 32 pixel window was divided. The block sizes corresponded to 1 × 1, 2 × 2, 4 × 4, and 8 × 8 pixels. In addition to these four standard deviation values, the standard deviation of these four values was computed. For the segments, these were calculated as averages of the area covered [18].
• 3. Textural features based on co-occurrence matrices of pixel values [19,20] extracted for grid elements and derived for segments as average values of the grid elements within segments:

  • $\begin{array}{cc}\text{Angular second moment}\hfill & {\displaystyle {\sum }_q{\displaystyle {\sum }_r{p}^2(q,r)}}\hfill \end{array}$
  • $\begin{array}{cc}\mathrm{Contrast}\hfill & {\displaystyle {\sum }_q{\displaystyle {\sum }_r{(q-r)}^2*p(q,r)}}\hfill \end{array}$
  • $\begin{array}{cc}\mathrm{Correlation}\hfill & {\displaystyle {\sum }_q{\displaystyle {\sum }_r(q*r*p(q,r)-{\mu }_x*{\mu }_y)}/}({\sigma }_x*{\sigma }_y)\hfill \end{array}$
  • $\begin{array}{cc}\mathrm{Entropy}\hfill & -{\displaystyle {\sum }_q{\displaystyle {\sum }_{r}p(q,r)*\log (p(q,r))}}\hfill \end{array}$
  • $\begin{array}{cc}\text{Local homogeneity}\hfill & {\displaystyle {\sum }_q{\displaystyle {\sum }_{r}p(q,r)/(1+{(q-r)}^2)}}\hfill \end{array}$

where:

• $p(q,r)=M(q,r)/Nt$
• M(q,r) = the co-occurrence matrix of the requantified pixel values q and r
• Nt = the total number of possible pairs in the image window
• μ_x ,σ_x = the mean and standard deviation of the row sums of the co-occurrence matrix
• μ_y ,σ_y = the mean and standard deviation of the column sums of the co-occurrence matrix
• The textural features based on co-occurrence matrices of pixel values were extracted in 4 directions in the extraction window: horizontally (0° angle), vertically (90°) and diagonally (45° and 135°). Pixel lag of 3 meters was applied in extracting these features on the basis of earlier study [7].

In addition, the following features were extracted from the ALS height data only:

• 4. Height statistics for the first and last pulses of all ALS points inside the field plot area or the segment area. These included mean, standard deviation, maximum, coefficient of variation, heights where certain percentages of points (5, 10, 20, ..., 95) had accumulated, and percentages of points accumulated at certain relative heights (5, 10, 20, ..., 95). Only points over 2 m in height were considered in computation of these variables. Finally, the percentage of points over 2 m in height was included as a variable.

For the estimation of forest attributes, all aerial image and ALS features were standardized to a mean of 0 and a standard deviation of 1. This was done because the original features had very diverse scales of variation. Without standardization, variables with wide variation would have had greater weight in the estimation, regardless of their correlation with the estimated forest attributes.

2.5. Selection of Features and Estimation of Forest Attributes

The k-nearest neighbor (k-nn) method was used for estimating the forest variables (e.g., [21,22,23]). The estimated variables were total volume of growing stock; the volume of Scots pine, of Norway spruce, and of deciduous species; basal area; mean diameter; and mean height. The value of k was set to 5 in both study areas, which was a compromise between the estimation accuracy and the averaging allowed in the estimation results. The k-nn estimation method typically has an increasing trend in accuracy when one raises the value of k from 1 to 10 (e.g., [23,24]), but large values of k typically result in retaining less of the original variation in the estimation results, as well as disappearance of the rare strata in the study material.

Euclidean distances were used to measure the closeness in the feature space, and the nearest neighbors were weighted with the inverse squared distances. The accuracy of the estimates was calculated via leave-one-out cross-validation by comparing the estimated forest variable values with the measured values (ground truth) of the field plots. The accuracy of the estimates was measured in terms of the relative root mean square error (RMSE) (see Equation (1)).

$$RMSE\%=100*\frac{RMSE}{\overline{y}}$$

(1)

where:

• $RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{({\widehat{y}}_i-y_i)}^2}}{n-1}}$
• y_i = measured value of variable y on plot i
• ŷ_i = estimated value of variable y on plot i
• $\overline{y}$ = mean of the observed values
• n = number of plots.

Automatic feature selection was carried out by means of a simple genetic algorithm presented by Goldberg [25] and implemented in the GAlib C++ library [26]. The reason for selecting this method was its success in an earlier study by Haapanen and Tuominen [27]. The GA process starts by generating an initial population of strings (chromosomes or genomes), which consist of separate features (genes). The strings evolve during a user-defined number of iterations (generations). This evolution includes the following operations: selecting strings for mating by applying a user-defined objective criterion (the better, the more copies in the mating pool), allowing the strings in the mating pool swap parts (cross over), causing random noise (mutations) in the offspring, and passing the resulting strings to the next generation.

In the present study, the starting population consisted of 300 random feature combinations (genomes). The length of the genomes corresponded to the total number of features in each step, and the genomes contained a 0 or 1 at position i, denoting the absence or presence of image feature i. The number of generations was 30. The objective variable to be minimized during the process was a weighted combination of relative RMSEs of k-nn estimates for mean volume, volume of Scots pine, Norway spruce volume, volume of deciduous species, mean diameter, and mean height, with mean volume having a weight of 50% and the remaining variables weighted at 10% each. Genomes that were selected for mating swapped parts with each other with a probability of 80%, producing offspring. Occasional mutations (flipping 0 to 1 or vice versa) were added to the offspring (with a probability of 1%). The strings were then passed to the next generation. The overall best genome of the current iteration was always passed to the next generation, also. Four successive steps (all including 30 generations) were taken, to reduce the number of features to a reasonable minimum. Only features belonging to the best genome in each step were included in the next step. Feature selection was run separately for both areas and each feature extraction unit (Grid, Seg350, and Seg1000). The list of selected features for all extraction units is presented in the appendix.

3. Results

In both study areas, the features extracted from square grid elements worked better in estimation of the forest attributes than the features extracted from image segments did (see Table 2 and Table 3). Furthermore, features from image segments derived with a minimum size of 350 m² performed better in the estimation than did features extracted from larger segments (minimum size: 0.1 ha).

Study area 2 had generally better estimation accuracy in comparison to the data sets of Study area 1. The main reason for this is probably the higher number of sample plots in Study area 2, which yields a higher number of potential nearest neighbors for each sample plot in the k-nn estimation.

There were large differences between the study areas in the estimation accuracy for volumes for the tree species groups. Typically, the highest estimation accuracy was seen with the dominant tree species. However, the volume of deciduous trees showed better estimation accuracy in comparison with the minority coniferous tree species group, since the presence of the deciduous trees is more easily recognizable in the aerial images.

Table 2. Estimation results for the feature sets: RMSE (%) relative to average (RMSE avg.) and standard deviation (RMSE std.) of Study area 1.

	Grid		Seg350		Seg1000
	RMSE Avg.	RMSE Std.	RMSE Avg.	RMSE Std.	RMSE Avg.	RMSE Std.
Total volume	27.8	43.0	34.0	52.6	36.6	56.7
Volume of Scots pine	74.2	59.6	77.1	61.9	99.9	80.4
Volume of Norway spruce	83.9	56.3	87.5	58.7	103.3	69.2
Volume of deciduous species	85.3	68.8	88.7	71.6	93.9	76.3
Basal area	25.8	49.8	30.1	58.1	29.8	57.7
Height	18.5	46.9	22.4	56.7	25.5	64.7
Diameter	25.5	57.2	27.7	62.1	32.0	71.9

Table 2. Estimation results for the feature sets: RMSE (%) relative to average (RMSE avg.) and standard deviation (RMSE std.) of Study area 1.

	Grid		Seg350		Seg1000
	RMSE Avg.	RMSE Std.	RMSE Avg.	RMSE Std.	RMSE Avg.	RMSE Std.
Total volume	27.8	43.0	34.0	52.6	36.6	56.7
Volume of Scots pine	74.2	59.6	77.1	61.9	99.9	80.4
Volume of Norway spruce	83.9	56.3	87.5	58.7	103.3	69.2
Volume of deciduous species	85.3	68.8	88.7	71.6	93.9	76.3
Basal area	25.8	49.8	30.1	58.1	29.8	57.7
Height	18.5	46.9	22.4	56.7	25.5	64.7
Diameter	25.5	57.2	27.7	62.1	32.0	71.9

Table 3. Estimation results for the feature sets: RMSE (%) relative to average (RMSE avg.) and standard deviation (RMSE std.) of Study area 2.

	Grid		Seg350		Seg1000
	RMSE avg	RMSE std	RMSE avg	RMSE std	RMSE avg	RMSE std
Total volume	29.6	43.1	32.9	48.1	34.8	50.6
Volume of Scots pine	125.2	75.7	138.5	82.3	137.0	81.6
Volume of Norway spruce	59.0	47.4	61.5	50.0	65.8	53.5
Volume of deciduous species	99.2	63.1	113.5	73.0	111.3	70.4
Basal area	25.3	50.5	26.9	53.6	28.0	55.7
Height	12.5	31.6	13.9	35.2	16.5	41.6
Diameter	19.8	41.2	22.7	47.3	25.2	52.3

Table 3. Estimation results for the feature sets: RMSE (%) relative to average (RMSE avg.) and standard deviation (RMSE std.) of Study area 2.

	Grid		Seg350		Seg1000
	RMSE avg	RMSE std	RMSE avg	RMSE std	RMSE avg	RMSE std
Total volume	29.6	43.1	32.9	48.1	34.8	50.6
Volume of Scots pine	125.2	75.7	138.5	82.3	137.0	81.6
Volume of Norway spruce	59.0	47.4	61.5	50.0	65.8	53.5
Volume of deciduous species	99.2	63.1	113.5	73.0	111.3	70.4
Basal area	25.3	50.5	26.9	53.6	28.0	55.7
Height	12.5	31.6	13.9	35.2	16.5	41.6
Diameter	19.8	41.2	22.7	47.3	25.2	52.3

There were 19 (Study area 1) or 14 (Study area 2) features selected for the final Grid sets, 19 or 17 for the Seg350 sets, and 12 or 18 for the Seg1000 sets. Of the selected features, the majority (61–79%) were based on the ALS data (height or intensity). In the Norway-spruce-dominated Study area 2, the proportion of features based on ALS height was larger, whereas in Study area 1 with its more even distribution of tree species, features based on ALS intensities were needed more often, probably for distinguishing among species.

Despite the larger number of aerial-photograph-based features available in Study area 2 (96 vs. 72 in Study area 1), their proportion in the final feature sets was smaller in the Grid case (21% vs. 37% in Study area 1) and only slightly greater in the Seg350 case (35% vs. 32%). For Seg1000, however, the proportion of aerial-photograph-based features was 39% in Study area 2 (in contrast to 33% in Study area 1). The total number of features selected also increased in moving from Grid level to Seg1000 level in Study area 2. In Study area 1, there was no trend in the number or properties of the selected features when one moved from Grid to Seg1000 level.

The need for study-area-specific feature selection can be seen from the above results. However, we wanted to see the effect of using features optimized for another geographical area and run the leave-one-out cross-validations for both study areas at the Grid level while using features optimized for the other area. This was possible since features based on the blue channel of aerial photographs, available for Study area 2 only, were not selected for the final Grid-based feature set of that area. Results are shown in Table 4. The results are slightly worse than those in Table 2 and Table 3, the only exception being deciduous species in Study area 1.

Another test was carried out at the Grid level, where the objective variable to be minimized during the GA process was the RMSE of mean volume, instead of the combined variable of the base case. The relative mean volume RMSEs were lowered by about three percentage points. The species-specific volumes of Study area 1 were heavily affected, with relative RMSEs rising by 21, 25, and eight percentage points for Scots pine, Norway spruce, and deciduous species, respectively. In the Norway-spruce-dominated Study area 2, the increases were 10, 2 and 14 percentage units.

Table 4. Estimation results (relative RMSE, %) with feature sets optimized for another study area, Grid-level results.

	Study area 1	Study area 2
Total volume	30.6	31.6
Volume of Scots pine	80.3	135.9
Volume of Norway spruce	87.5	64.8
Volume of deciduous species	80.0	113.4
Basal area	26.7	26.8
Height	18.6	13.1
Diameter	25.2	20.2

Table 4. Estimation results (relative RMSE, %) with feature sets optimized for another study area, Grid-level results.

	Study area 1	Study area 2
Total volume	30.6	31.6
Volume of Scots pine	80.3	135.9
Volume of Norway spruce	87.5	64.8
Volume of deciduous species	80.0	113.4
Basal area	26.7	26.8
Height	18.6	13.1
Diameter	25.2	20.2

We took a closer look at the errors in different volume classes, in order to gain fuller understanding of the effect of forest structure on the results. Figure 3 shows volume estimation errors by volume class at the Grid level. It seems that the greater volume estimation errors in Study area 2 are due to the field plots with the highest volumes, obviously lacking suitable neighbors. The last two volume classes were completely absent from Study area 1. Lack of suitable neighbors was seen also in the results obtained with feature sets optimized for estimation of mean volume only: the mean volume RMSEs decreased, but those of species-specific volumes rose, meaning that there were more eligible neighbors with similar volumes but the actual forest composition differed.

Figure 3. Volume estimation errors by volume class at Grid level.

Figure 3. Volume estimation errors by volume class at Grid level.

4. Discussion and Conclusions

In this study, the general level of estimation accuracy at the plot level was similar to, or slightly lower than in studies utilizing ALS and aerial image data in similar conditions (e.g., [3,28,29]). The slight differences in estimation accuracy are, presumably, caused mainly by differences in the properties of the study areas, such as variation in stand structures and tree species compositions. Despite the lower pulse density of ALS data in Study area 2, the relative basal area, height, and diameter errors were smaller, and the relative mean volume error was only slightly greater than in Study area 1.

Despite the theoretical advantages of the segment-based approach, the features extracted for segments did not perform well in the estimation procedure. The main reason can be assumed to be the study setting: the field data were measured per sample plot and not per segment produced. If the automatic segmentation produces segments that are homogeneous enough, this should not be a significant problem. However, our results suggest that the segmentation does not seem to delineate segments optimally from the standpoint of feature extraction for stand variable estimation. Furthermore, automatic segmentation often produces segments that are irregularly shaped—i.e., not compact. In forest stands with large trees, the segment borders are typically found in gaps between trees, in which case the variation within the segments may be more significant than the variation between segments. On the other hand, using geographically larger segments typically resulted in lower estimation accuracy in comparison to other feature sets, which indicates that the larger the units are the more internal variation they have. Thus, increasing the size of feature extraction unit seems to have a negative effect on the estimation accuracy. In a study by Pekkarinen and Tuominen [11] that utilized only aerial photograph data, the results were similar in the size of the segments, although the results were inconclusive in relation to the superiority of the plot- or segment-based methods. In that study, the segment-based approach seemed to retain better the original spectral variation of the aerial photographs, as there was less averaging in the segment-based image features [11].

From another viewpoint, the comparison carried out in this study is one way to evaluate the success of automatic stand delineation. With auxiliary data having a high correlation with the field data, the estimation errors should be close to each other regardless of the extraction unit.

The feature selection process was carried out separately for the two study areas and each extraction unit set. Although the randomness inherent in the applied genetic algorithm method results in somewhat random final feature sets, there were logical differences between the study areas, with ALS-height-based features playing a greater role in the area with one clearly dominant tree species (Study area 2, with Norway spruce), whereas ALS-intensity-based features aided in discriminating between tree species in the more heterogeneous area. This result is backed up by the fact that there were a larger number of aerial-photograph-based features available in the Norway-spruce-dominated area, yet their proportion was smaller or about the same as in the area with even proportions of Scots pine, Norway spruce, and deciduous tree species. Furthermore, using features optimized for another area resulted in slightly larger errors than obtained with features optimized for the area in question. The tree-species-specific estimates for area 2, with its domination by Norway spruce, were affected most.

The estimation accuracy always includes a component of error due to the errors in field measurements and in modeling of the variables—such as volume—on the basis of the individual tree measurements. The errors in the field measurements and the tree models used were not taken into account in assessment of the estimation accuracy; the field reference data were considered an exact representation of local forest attributes.

From the results of this study, the most feasible inventory procedure utilizing ALS and aerial image data seems to be the following: (1) extraction of ALS- and aerial-photograph-based features for the systematic grid elements; (2) area-specific optimization of the feature combination to be employed; (3) estimation of forest variables for the systematic grid elements; (4) automatic segmentation utilizing ALS height, ALS intensity, and aerial photographs; (5) derivation of the estimates for image segments on the basis of the estimates of grid elements; and (6) manual combination of image segments for purposes of deriving spatial units for utilisation in forest management.