Comparison of grid-based and segment-based estimation of forest attributes using airborne laser scanning and digital aerial imagery

Remote Sens. 2011, 3, 945-961; doi:10.3390/rs3050945 Remote Sensing ISSN 2072-4292 www.mdpi.com/journal/remotesensing Article Comparison of Grid-Based and Segment-Based Estimation of Forest Attributes Using Airborne Laser Scanning and Digital Aerial Imagery Sakari Tuominen 1,* and Reija Haapanen 2 1 Finnish Forest Research Institute, Metsäntutkimuslaitos, PL 18, 01301 Vantaa, Finland 2 Haapanen Forest Consulting, Kärjenkoskentie 38, 64810 Vanhakylä, Finland; E-Mail: [email protected] * Author to whom correspondence should be addressed; E-Mail: saka[email protected]; Tel.: +358-10-211-2167; Fax: +358-10-211-2202. Received: 2 February 2011; in revised form: 4 March 2011 / Accepted: 2 May 2011 / Published: 12 May 2011 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 OPEN ACCESS Remote Sens. 2011, 3 946 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. Keywords: forest inventory; airborne laser scanning; aerial photography; image segmentation 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–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/m2) 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 Remote Sens. 2011, 3 947 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. Remote Sens. 2011, 3 948 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 Remote Sens. 2011, 3 949 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 m2 and 0.1 ha, were carried out in both study areas. Remote Sens. 2011, 3 950 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. (a) (b) (c) 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 Remote Sens. 2011, 3 951 whose minimum size was set as 350 m2. 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:  Angular second moment   q r rqp ),( 2  Contrast    q r rqprq ),(*)( 2  Correlation )*(/)*),(**( yx q r yx rqprq      Entropy   q r rqprqp )),(log(*),(  Local homogeneity    q r rqrqp ))(1/(),( 2 where: NtrqMrqp /),(),(  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. Remote Sens. 2011, 3 952 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–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)). y RMSE RMSE *100% (1) where: 1 ) ˆ ( 1 2      n yy RMSE n iii yi = measured value of variable y on plot i ŷi = estimated value of variable y on plot i 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 Remote Sens. 2011, 3 953 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 Tables 2 and 3). Furthermore, features from image segments derived with a minimum size of 350 m2 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 Remote Sens. 2011, 3 960  Combined standard deviation of pixel blocks (1 × 1, 2 × 2, 4 × 4, 8 × 8) of aerial image NIR band  Angular second moment (45° angle) of aerial image red band  Angular second moment (45° angle) of aerial image green band  Homogeneity (135° angle) of aerial image green band  Maximum of first pulse hits  Standard deviation of first pulse hits (below 2 m hits excluded)  Height, where 20% of first pulse hits have been accumulated (below 2 m hits excluded)  Height, where 80% of first pulse hits have been accumulated (below 2 m hits excluded)  Percentage of first pulse hits below 70% of maximum height (below 2 m hits excluded)  Percentage of last pulse hits above 2 m height  Height, where 20% of last pulse hits have been accumulated (below 2 m hits excluded)  Percentage of last pulse hits below 95% of maximum height (below 2 m hits excluded) Seg1000  Average of ALS height  Standard deviation of ALS intensity  Standard deviation of 2 × 2 pixel blocks of ALS intensity  Contrast (135° angle) of ALS height  Standard deviation of aerial image NIR band  Contrast (135° angle) of aerial image NIR band  Contrast (90° angle) of aerial image red band  Contrast (135° angle) of aerial image green band  Height, where 90% of first pulse hits have been accumulated (below 2 m hits excluded)  Height, where 10% of last pulse hits have been accumulated (below 2 m hits excluded)  Percentage of last pulse hits below 30% of maximum height (below 2 m hits excluded) Study area 2. Grid  Average of ALS height  Average of ALS intensity  Contrast (135° angle) of ALS height  Entropy (0° angle) of ALS height  Average of aerial image NIR band  Entropy (135° angle) of aerial image green band  Entropy (90° angle) of aerial image NIR band  Height, where 10% of first pulse hits have been accumulated (below 2 m hits excluded)  Height, where 40% of first pulse hits have been accumulated (below 2 m hits excluded)  Height, where 90% of first pulse hits have been accumulated (below 2 m hits excluded)  Percentage of first pulse hits below 80% of maximum height (below 2 m hits excluded)  Percentage of last pulse hits above 2 m height Remote Sens. 2011, 3 961 Seg350  Average of ALS height  Angular second moment (135° angle) of ALS intensity  Homogeneity (90° angle) of ALS height  Standard deviation of 4 × 4 pixel blocks of aerial image blue band  Average of aerial image NIR band  Angular second moment (90° angle) of aerial image red band  angular second moment (0° angle) of aerial image NIR band  Entropy (0° angle) of aerial image green band  Entropy (0° angle) of aerial image NIR band  Percentage of first pulse hits above 2 m height  Std of first pulse hits (below 2 m hits excluded)  Height, where 10% of first pulse hits have been accumulated (below 2 m hits excluded)  Height, where 40% of first pulse hits have been accumulated (below 2 m hits excluded)  Height, where 80% of first pulse hits have been accumulated (below 2 m hits excluded)  Height, where 30% of last pulse hits have been accumulated (below 2 m hits excluded)  Percentage of last pulse hits below 10% of maximum height (below 2 m hits excluded)  Percentage of last pulse hits below 90% of maximum height (below 2 m hits excluded) Seg1000  Average of ALS height  Average of ALS intensity  Angular second moment (0° angle) of ALS intensity  Entropy (135° angle) of ALS height  Homogeneity (0° angle) of ALS height  Homogeneity (90° angle) of ALS height  Average of aerial image red band  Average of aerial image NIR band  Standard deviation of aerial image NIR band  Average of aerial image blue band  Standard deviation of 8 × 8 pixel blocks of aerial image blue band  Homogeneity (0° angle) of aerial image green band  Percentage of first pulse hits above 2 m height  Height, where 50% of first pulse hits have been accumulated (below 2 m hits excluded)  Height, where 80% of last pulse hits have been accumulated (below 2 m hits excluded)  Percentage of last pulse hits below 10% of maximum height (below 2 m hits excluded)  Percentage of last pulse hits below 20% of maximum height (below 2 m hits excluded) © 2011 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/).