A Two-Level Approach for Species Identification of Coniferous Trees in Central Ontario Forests Based on Multispectral Images

IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING, VOL. 8, NO. 4, APRIL 2015

1487

A Two-Level Approach for Species Identification of Coniferous Trees in Central Ontario Forests Based on Multispectral Images Jili Li, Baoxin Hu, Member, IEEE, and Murray Woods Abstract—This study aims to provide detailed spatial information of valuable tree species to support improved management of winter habitat of white-tailed deer. To achieve this, we proposed a novel approach using information from two spatial scales and a suite of methods for analysis and classification of remotely sensed data. High-spatial resolution, multispectral images were employed to test the proposed method. A new structure-based remote sensing feature [local binary pattern (LBP) index] was developed and proved to be effective for species classification. A simple but effective fusion approach based on information entropy theory was proposed to incorporate features derived from different methods and their uncertainties. Based on tenfold cross validation, an overall accuracy (OA) of 77% was obtained for the classification of three tree species groups. The proposed approach has high potential to improve species mapping for operational ecological modeling. Index Terms—Entropy, forestry, image processing.

I. I NTRODUCTION

T

REE SPECIES is arguably the most important and widely sought piece of information in forest inventory and management, since it serves as the basis for estimating a variety of attributes of individual trees or tree stands, such as volume, biomass [1], and fiber quality. Information on the distribution of certain species is also needed for wildlife habitat mapping and for conserving the quality of forest ecosystems. For this study, our objective was to provide information of valuable species of individual trees to facilitate the management of the winter habitat of white-tailed deer (Odocoileus virginianus). Although abundant food sources make almost any forested area suitable for white-tailed deer during summers, the deer concentrate in areas with certain conifer species that provide food and shelter from cold, storms, and deep snow during winters (commonly called deer yards). It is important to know the distribution of species in forests to accurately assess the winter habitat suitability. Manuscript received November 03, 2014; revised March 27, 2015; accepted April 09, 2015. Date of publication April 29, 2015; date of current version May 26, 2015. This work was supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada. J. Li and B. Hu are with the Department of Earth and Space Science and Engineering, York University, Toronto, ON M3J1P3, Canada (e-mail: [email protected]; [email protected]). M. Woods is with the Ontario Ministry of Natural Resources and Forestry, North Bay, ON P6A 2E5, Canada (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/JSTARS.2015.2423272

Despite the long history of research and development in tree species identification from remotely sensed data, automatic species classification remains challenging and an ongoing active research area. Traditionally, medium and low spatial resolution satellite imagery is used to classify forest cover types mostly based on their spectral information (e.g., [2]). To classify individual species, researchers have exploited hyperspectral imagery to take advantage of the subtle difference among species in their reflectance spectra (e.g., [3]–[5]). However, the spatial resolutions of most of the available hyperspectral images are not very high (greater than 1 m) and thus most of the classification methods are pixel-based. Even though progress has been made in species classification using hyperspectral imagery, the main applications are still focused on stand or tree-group level. Identification of individual species at individual tree levels is made available with the development of high-spatial resolution imagery and small footprint light detection and ranging (LiDAR) data [6]. Species identification of individual trees has been conducted using object-based approaches [7], which have been demonstrated to be effective in analyzing remotely sensed imagery with high spatial resolution [8]. Compared with the traditional pixel-based approaches, using an object-based method, tree crowns are first detected, segmented, or delineated from high spatial resolution remotely sensed data; various features are then extracted to characterize each tree object and used for classification. Even though the success of species identification with airborne LiDAR data is widely accepted, the data acquisition, processing, and management of LiDAR data are far more expensive than those of imaging data. Studies have been conducted to exploit airborne (passive) optical imagery in characterizing the structural information of tree canopies (e.g., [9] and [10]) as a cost-effective alternative to LiDAR data. In the literature, remote sensing features used within the context of species classification using high spatial resolution imagery are mostly spectral or statistical texture measures. Leckie et al. [11] presented a species classification study in oldgrowth conifer stands of western Canadian forests using high spatial resolution imagery. In their study, the mean reflectance of the sunlit part of a tree crown was adopted as a unique spectral feature to classify the species of segmented tree crowns. Although simple spectral classification with carefully selected species can produce positive results, they noted tremendous variability within crowns of the same species in old-growth stands and overlap in spectral characteristics among crowns of different species, which made species identification difficult.

1939-1404 © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

1488

IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING, VOL. 8, NO. 4, APRIL 2015

Brandtberg [12] exploited different remote sensing features derived from optical imagery for individual species classification. Their features included the spectral values, concaved and curved contours, and radial pattern of tree crowns. An OA of 64% was achieved to classify nine tree species. In addition, the authors indicated the difficulty in identifying the optical features that are systematically different among species on their study area. Waser et al. [13] present an approach employing two groups of explanatory variables to classify seven tree species in Switzerland from airborne digital scanner-40 (ADS40) (0.5 m spatial resolution) images. The developed spectral and geometric variables from the optical image yielded classification results of 70%–80% accuracies depending on their test sites. Immitzer et al. [14] used high spatial resolution WorldView-2 satellite data (0.5 m at the panchromatic band and 2 m at multispectral bands) to classify ten tree species in a temperate forest in Austria. The spectral variability among species was analyzed and the object-based classification employing the multispectral features was proved better than the pixel-based approach. An OA of 82% was achieved while the class-specific accuracies ranged from 37% to 94%. In the context of high spatial resolution optical imagery, tree structures may be described by the textural features of tree crowns [15]. Textural information has been widely used in remote sensing of vegetation. Most of the existing studies are focused at the stand or landscape scales where the spatial distribution of trees forms the dominant texture [16], [17]. At the crown level, in [18] and [19], the gray-level cooccurrence matrix (GLCM) and its corresponding statistical measures, such as energy, contract, entropy, and homogeneity, were utilized for the identification of multiple tree species based on imagery. The obtained classification accuracies varied between 0.4 and 0.87 for different species. Zhang and Hu [10] developed a longitudinal profile feature describing the side view of tree shapes. The new feature was combined with a decision tree classifier to improve urban tree species identification from high spatial resolution airborne imagery. However, there are few studies exploited the textural information to characterize the spatial distribution of leaves and branches within individual tree crowns (ITCs). The aforementioned methods have been mostly limited in the use of spectral information and statistical textural measures, which reveals there seems to be space left for improvements on the identification of tree species and in the development of new effective features and methodology. Different tree species tend to have different foliage structures in both vertical and horizontal directions [20]. The spatial distribution of foliage within a crown also account for major differences among species, such as foliage gaps and clumps/clusters. These physical structural features normally result in primitives of textures on high resolution remotely sensed imagery. The application of structural textural measures is useful to characterize the porosity and spatial distribution of gaps within a crown and has been proved successful in computer vision and pattern recognition [14]. Then, our first research question was whether vertical and horizontal structures of a tree crown can be described by structural texture measures which could be used to improve tree species classification.

Even though attributes of ITCs have been well investigated in species classification, the contextual information related to individual trees remain under studied. It is commonly known and accepted that, trees of a certain species tend to grow together on the landscape. Remotely sensed imagery at the crown-level resolution (e.g., 2–4 m) may be able to provide the contextual information and support individual tree classification. As an example, an individual deciduous tree maybe misclassified as one of coniferous species based on the information provided by high spatial resolution imagery, due to either the similarity in features among species or the anomaly of this tree in features compared with other trees of the same species. However, if it is known that this tree is in a stand dominated by deciduous species, the classification result may be different. In the literature, contextual information is often used in postprocessing of classification results based on the majority rules, but is not used intelligently. A tree that is in a stand dominated by deciduous species does not necessarily belong to a deciduous species group. The second research question we were to address in this study is how integrate the information at the level of tree stands in the crown-based classification. To address the first research question, we investigated a new textural descriptor named local binary pattern (LBP) and its rotation-invariant version. A new index feature based on the LBP was designed to describe the tree crown texture, and the effectiveness of the index was demonstrated by our experiment. The proposed texture measure is the core contribution and novelty of this paper. The second research question was conquered by introducing a two-level classification framework and linking the two level classification results through measures of classification uncertainty. The proposed new framework is also original and proved to be effective. II. S TUDY A REA AND DATA The study area is located in the district of Parry Sound, Ontario, Canada (45.33N, 80.03W). Forests in the study area are mostly composed of conifers, with a large population of eastern hemlock (Tsuga canadensis), pine, and spruce. Aerial orthophotos over the study area were obtained using a Leica ADS40 in the summer of 2009. These orthophotos have a spatial resolution of 0.4 m. A subset of the imagery is shown in Fig. 1. There are four spectral bands in the near-infrared, red, green, and blue regions, respectively. The delivered image data were processed to produce at-surface reflectance recorded in 16-bit integer values. A total of six forest sites were selected to test our methods. Within the spatial range of these six sites, 181 sample trees containing eight coniferous species were selected to establish a training/testing dataset for this study. The six sites were selected according to the forest accessibility, while each site was ensured to cover at least two dominant or codominant species in this study. Sample trees were all in good growth condition and had normal crown shapes and sizes. The coniferous species are eastern hemlock, eastern white pine (Pinus strobus), red pine (Pinus resinosa), black spruce (Picea mariana), balsam fir (Abies balsamea), european larch (Larix decidua), and eastern white cedar (Thuja occidentalis). Species of these sample

LI et al.: TWO-LEVEL APPROACH FOR SPECIES IDENTIFICATION OF CONIFEROUS TREES

1489

Fig. 1. Initial multispectral ADS40 image with 0.4 m spatial resolution. Fig. 2. Diagram shows the workflow of the proposed method.

trees were determined during a field inspection in May 2014. These species were grouped into three classes: 1) high value species (HS) including eastern hemlock and cedar; 2) medium value species (MS) including larch, white spruce, black spruce, and balsam fir; and 3) low value species (LS) including white pine and red pine. The rationale behind the grouping of these species is twofold. From the perspective of ecology, the species in each group have similar ecological values in term of the suitability for the winter habitat of white tailed deer. Ecologists are content with the separation of the species among groups for white-tailed deer habitat mapping (Dr. Naylor, personal communication). In term of image classification, for a given group, some of the tree species have similar spectral and structural features in the remote sensing data and thus it is hard to separate them. For example, larch, spruce, and fir are mostly belong to subcanopy layer trees and usually in narrow cone shape and pointing tree tops, while hemlock and cedar both have denser leaf and branches, and green tones. By grouping these species together, the within-class variations in term of the features used in this study were presumably small and the between-class variations were large. Even though the species in one group may be clustered together in the feature domain considering all of the features used in this study, for an individual feature, the species within a group may not necessarily be similar. The number of sample trees within each class was 91, 55, and 35, respectively. Hardwood trees in the study area were treated as one group as they are not essential in forming deer yards.

III. M ETHODS In this study, we propose a simple but effective two-level classification approach: pixel-wise classification at crown-level (4 m pixel) for broad land-cover classes: coniferous and deciduous canopies, bare earth, water bodies, and grass (level 1); object-based classification of ITCs from high resolution images (0.4 m pixel) for three types of coniferous trees: HS, MS, and LS for deer yards (level 2). Classification results from these two levels including their uncertainties were then considered for the final decision through a set of logic-based decision rules. In our method, the contextual information associated with individual trees was employed through the classification at the first level

and the detailed structural information within tree-crowns was exploited at the second level. The proposed method consisted of the following four stages (Fig. 2): pixel-level classification, ITC segmentation, feature extraction, and object-based classification and decision fusion.

A. Level 1: Pixel-Level Classification The high-spatial resolution image was first downsized by pixel aggregation to have a spatial resolution of 4 m × 4 m, so that the resulting pixels represent relatively homogeneous stands. The choice of downsizing to 4 m was made by considering the general size of typical tree groups/stands in our study area and also the possible grid size that may affect deer movement. It is worth mentioning that any remotely sensed imagery at a similar resolution is suitable. It was expected that a simple and fast classification method based on the spectral reflectance in the visible and near-infrared bands could be selected to identify broad forest cover types (i.e., coniferous and deciduous). We tested both the maximumlikelihood (ML) and support vector machine (SVM) methods for the first-stage classification. Both methods generated similar results, but as expected SVM was slightly better (about 3% improvement in the first stage classification validation). However, considering the large number of pixels that needed to be classified at the first stage, we decided to use the simple and fast ML method. In addition, we believe that ML could generate consistently good classification results considering the assumption that the features for each pixel class follow multivariant Gaussian distributions is generally reasonable. ML classification was used to classify individual pixels on the downsized image into five cover types: coniferous canopy, deciduous canopy, bare earth, water, and grass. Training areas for each class with homogenous pixels were visually interpreted by forest experts from the original high resolution imagery (0.4 m). We applied all of the four multispectral components as the input features for the classification process. The training pixels were from several region of interest across the image and the total numbers of training pixels were 8,656, 6183, 7512, 8,438, and 3385 for the above five cover types, respectively.

1490

IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING, VOL. 8, NO. 4, APRIL 2015

The classification accuracy was evaluated based on a number of 100 stratified randomly selected sample pixels per class (a total 500 test samples). The test samples were not overlapped with the training samples. The resultant classification image was resized back to the original spatial resolution (0.4 m) for later analysis. B. Level 2: ITC Segmentation In our previous research [21], we developed an effective automatic ITC segmentation algorithm by utilizing a multiscale filtering and segmentation (MFS) technique. The MFS algorithm was also used on the original 0.4 m resolution image (green band) in this research, with a few modifications on the masking procedure. We used the normalized difference vegetation index (NDVI) to separate vegetation and nonvegetation areas by a single threshold and the grass pixels classified in level 1 to separate grass and tree areas within the vegetation cover. The determined nontree area was masked out following similar morphological processing as described in [21]. In addition, we used a morphological closing operator to shrink each segment by two pixels along the boundary. It is worth mentioning that besides the reflectance of in green band, we tested other spectral features, such as reflectance in the near-infrared band and brightness values generated by the principal component analysis and there was not much difference in the segmentation results. After the ITC segmentation, each segment was recorded as an object representing an ITC. When overlaying the level 1 classification map over the ITC segmentation map, it was possible that a segment covered several pixels belonging to different classes. We defined a measure of class mixture S based on Shannon’s information entropy S(x) = −

K 

pk (x)log2 pk (x)

(1)

belonging to one of the five class types defined at level 1 or UL1 , we were able to further remove segments labeled as deciduous canopy and other nonconiferous classes. In other words, only segments with a label of “coniferous canopy” or “ UL1 ” were retained. We kept the segments of UL1 because we assumed that information from individual tree classification can be used to further determine species of the segments whose class could not be determined through the level 1 classification. C. Level 2: Feature Extraction At this stage, various remote sensing features were extracted for each object (i.e., ITC) and used as input to a selected classification algorithm. In this study, two categories of features, statistical and structural, were derived. The features in the first category include first-order and second-order statistics (based on GLCM) of the original reflectance values. In order to complement the statistical features and characterize detailed structures of a tree crown, a new feature based on the LBP method was designed. It was expected that the involvement of the structure-based LBP feature could improve the classification. The features and feature extraction approaches are described as follows. 1) Statistical Features: For each tree crown, the mean and standard deviation of its pixel values in each of the four spectral bands were calculated. These statistics derived from the original reflectance values revealed the most obvious spectral differences among tree species. The horizontal/vertical distribution of leaves and branches within a tree crown constitutes the tree’s texture. To characterize the image texture, the gray scale cooccurrence matrix (GLCM) was used. It was calculated with pixel shifts of a series of distances d [22]. Two statistical features, Homogeneity and angular second moment (ASM), were then calculated from GLCMs for each spectral band

k=1

where pk (x) is the probability of class k in the segment x, k represents the kth class, and K is the number of classes which equals to 5. pk can be estimated as the ratio of the number of pixels in class K to the total number of pixels in the segment. The entropy reflected the impurity of a given segment. It is clear from (1) that S reaches to the highest value when every class has an equal probability, and the lowest value when only one class present in the segment. Therefore, a large value of S implied that the segment/ITC consisted of a very mixed classes determined at level 1; a small value of S indicates that the segment/ITC was dominant by one class. If the former scenario occurred, it was hard to make a decision only based on the level 1 classification due to the class mixture and thus we defined the segment as an unknown segment at level 1. We separated the known and unknown segments by using a predefined threshold TL1 on S. If S(x) > TL1 , the segment x was regarded as an unknown object at level 1 and given a unique label of UL1 . Otherwise, we assigned x the unique class label corresponding to the class occupying the maximum proportion of pixels within the segment. For example, if one class has more than 50% coverage in the segment, it will be assigned to this segment. Once each segment was given a unique label either

Homogeneity =

N N   i=1 j=1

ASM =

N  N 

1 1 + (i − j) p(i, j)

2 p(i, j)

2

(2)

(3)

i=1 j=1

p(i, j) =

GLCM (i, j) N  N  GLCM (i, j)

(4)

i=1 j=1

where N is the quantization level of the GLCMs, and i and j represent the row and column numbers in the GLCMs, respectively. In this study, N was set to 16 and d was set to 3. 2) Structural Features: The above-mentioned textural measures described the statistical relationships between a pair of pixels within a given crown, which reflected relationships among tree elements. However, the spatial distribution of these tree elements was not sufficiently described. Specifically, these statistical textual measures may characterize the gap fractions within a crown, but not the spatial distributions of the gaps. In order to characterize the detailed spatial distributions of tree elements (e.g., the contexture relation between withincrown gaps and neighboring foliage; the directional distribution

LI et al.: TWO-LEVEL APPROACH FOR SPECIES IDENTIFICATION OF CONIFEROUS TREES

1491

Fig. 4. Nine uniform patterns and one nonuniform patter that can occur in an image.

Fig. 3. Calculation and mapping of rotation invariant LBP. (a) The thresholded neighborhood. (b) Mapped binary value of the thresholded neighborhood. (c) An arbitrary number of binary shifts. (d) Matched one of the 36 local binary patterns.

of the clustered foliage of conifers) or the micro-primitives of textures (e.g., edges, spots, curved corners, and flat local areas), and thus improve species classification, we exploited a structure-based LBP method. A detailed description of this method can be found in [23]. The LBP method is briefly reviewed in the following paragraphs. LBP is a mathematical operator first introduced [24] as a shift invariant measure for local image texture. The gray value of the center pixel in a sliding window is considered as a threshold to code surrounding neighborhood pixels. Given a central pixel with a gray value of gc , a binary pattern number (a set of 0 or 1) and its LBP code are calculated by comparing gc with its eight-neighbor pixel values using LBP =

8 

2p f (gp − gc )

(5)

p=1

 f (gp − gc ) =

1, 0,

if gp − gc ≥ 0 if gp − gc < 0

(6)

where f is a function of gp − gc , p = 1, 2, . . . , 8, and gp is the pixel value of the pth neighbor. The bordering pixels of an object were not considered and count for the calculation of LBP. The LBP operator produces 28 different output values, corresponding to the 28 different binary patterns that can be formed by the eight neighbor pixels. When the image is rotated, the gray value gp correspondingly moves around gc . Rotating a particular binary pattern naturally results in a different LBP value. To remove the rotation effect, binary values of the thresholded neighborhood are mapped into an 8-bit series in clockwise order [Fig. 3(b)]. An arbitrary number of binary shifts is then performed [Fig. 3(c)], until the series matches one of the 36 different patterns [Fig. 3(d)] of “0” and “1” that 8-bit series can form. The index of the matching pattern is then used as a rotation invariant LBP defined as LBPri . A total of 36 unique rotation-invariant LBPs can occur in an image. It has been found in [24] that the occurrence frequencies of the 36 individual patterns varied greatly, and this was also the case in our study case. Furthermore, a variable U is defined as the number of spatial transitions between 0 and 1 in the LBPri U (LBPri ) = |f (g8 − gc ) − f (g1 − gc )| +

8 

|f (gp − gc )

p=2

− f (gp−1 − gc )|.

(7)

A pattern that has a U value equal to or less than 2 is defined as uniform, and the others are defined as nonuniform. For example, the pattern 11100111 (two transitions) is uniform, whereas the pattern 11001001 (four transitions) is nonuniform. Patterns are called “uniform” because of their uniform neighbor structure which contains very few spatial transitions. Ojala et al. [24] defined nine uniform patterns (No. 1–9, Fig. 4) and 27 nonuniform patterns that can occur in an image, and the nonuniform patterns were grouped into a single pattern, giving an extra code No. 10. These uniform and nonuniform patterns quantify the occurrence statistics of individual rotation invariant patterns corresponding to certain micro-primitives of textures in an image, therefore, these patterns can be considered as structural feature detectors. Different crown types may be discriminated through these patterns. For example, some crowns have more No. 1 patterns with bright spots, which indicate that tree elements are surrounded by gaps; others may have more No. 4–6 patterns with curved edges or corners, which describe the transaction line between branch boundaries and within crown gaps. Furthermore, a uniform rotation-invariant pattern LBPriu2 can be formally defined as  LBPri , if U (LBPri ) ≤ 2 riu2 LBP = (8) 10, otherwise. Superscript riu2 reflects the use of rotation-invariant uniform patterns with U values of at most 2. The mapping from LBPri to LBPriu2 , which has ten distinct output values, can be implemented with a look-up table with 28 elements [24]. From LBPriu2 , nine rotation-invariant uniform and one grouped nonuniform patterns are derived. For an individual tree, the percentages of pixels of these ten patterns can be expressed by a normalized histogram with ten bins, and the sum of the frequencies of all patterns is equal to 1. After the calculation of LBPriu2 for each pixel (i, j) in a given tree crown, the tree can be described by a histogram generated via Hist(k) =

N  M 

F (LBP riu2 (i, j), k),

i=1 j=1

F (LBP

riu2

(i, j), k) =



1, 0,

k ∈ [1, 8]

(9)

ifLBPriu2 (i, j) = k (10) otherwise

where F is a function of LBPriu2 and k, N is the number of rows of the image, and M is the number of columns. Each bin in the histogram represents a LBP. Fig. 5 presents the mean frequencies of uniform (No. 1–9) and nonuniform (No. 10) LBPs of the three groups of tree species (i.e., HS, MS, and LS) to be classified in this study. As an example, the LBPs were calculated using all training samples in the green band.

1492

IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING, VOL. 8, NO. 4, APRIL 2015

Fig. 5. Averaged LBP histogram of three coniferous species groups. HS, high value species; LS, low value species; MS, medium value species.

A quick visual inspection of Fig. 5 suggests that patterns No. 4, No. 5, No. 9, and No. 10 have better separability among species than the other patterns. A new index was then designed to reveal and amplify the potential difference of these patterns among species LBP I =

LBP riu2 (u5) − LBP riu2 (u9) LBP riu2 (u5) + LBP riu2 (u9)

Fig. 6. Crown images of a hemlock (left) and a balsam fir (right) tree, showing different gap distribution and within-crown textural difference with LBPI values of 0.06 for the hemlock and 0.82 for the balsam fir.

(11)

where u5 and u9 represent the uniform patterns No. 5 and No. 9 in the LBPriu2 output histogram. The use of local binary index feature (LBPI) reduced the number of original LBP features from 10 to 1, while preserving the most discriminative power for classifying tree species. For each spectral band, LBPI was calculated. According to the definition (11), a large LBPI value of a tree crown indicates the crown has more smoothly increased brightness regions and linear textures (pattern #5), more homogeneous areas, or less dark spots (pattern #9), and vice versa. For example, several dark spots (i.e., surrounding eight-neighbors have higher green values than the center pixel) and few linear features are visible on a hemlock image crown in Fig. 6. In contrast, the brightness values of the balsam fir image crown is smoothly increasing from border to the treetop with less dark spots but more smoothed line features than the hemlock crown. The textural difference can be distinguished by their LBPI values: 0.06 for the hemlock and 0.82 for the balsam fir. In addition, the mean and standard deviation of LBPI for the three species groups are displayed in Fig. 7. The location of the means for the three classes is distinguishable, but considering the standard deviations, there are overlaps in the LBPI values among these classes. To determine if the difference in the mean values among these classes were significant enough for their separation, we carried out unpaired two-sample t-test. The null hypothesis was that each of the paired two species groups (e.g., HS and MS) had equal means, and the significant level of 5% was used. The null hypothesizes were all rejected for the pairs of {HS, MS}, {HS, LS}, and {MS, LS} species groups, indicating the population means were not the same at a significant level of 5%. D. Level 2: Classification To measure classification uncertainty, it is required that the selected classification algorithm has an embedded mechanism to provide a posterior classification or membership probability.

Fig. 7. Mean (dot) and standard deviation (error bar) values of LBPIs calculated from the 181 sample trees. HS, high value species; MS, medium value species; LS, low value species.

There are several classification algorithms having the capability such as ML, SVM, and fuzzy logic classification (a softclassification method). In this research, we employed the SVM which has been successfully used in various pattern recognition studies and often performs better than other classification algorithms [25]. The basic idea of SVM is to map a given set of binary labeled training data to a high dimensional feature space and determine an optimal decision hyperplane with the largest margin to separate training samples. The output of a SVM classification can be a vector that contains class labels for each object, or a probability matrix that contains probability estimates for each object to belong to each class. The classification probability is calculated from the decision values by an improved implementation [26] of Platt’s a posteriori probabilities [27], [28]. The class probabilities can be extended to the multiclass case by combining all binary classifiers class probability outputs [29]. Detailed description of the SVM algorithm can be easily found in literature (e.g., [30]). The Gaussian radial basis function was selected in this research, as it is commonly suggested in most of the other remote sensing classification studies. The accuracy of SVM is affected by the selection of two parameters: C and r. In this research, the grid search approach was applied. It searches for every parameter pair among a predefined range and with a predefined step interval. The parameter pair producing the best classification accuracy (by cross validation) is chosen. A multiclassification problem in SVM is normally decomposed into a few binary classification problems. In this paper, we used the

LI et al.: TWO-LEVEL APPROACH FOR SPECIES IDENTIFICATION OF CONIFEROUS TREES

1493

TABLE I D ECISION F USION FROM T WO -S CALE C LASSIFICATION L ABELS

C, coniferous canopy; unknown C, unknown conifers; HS, high value species; MS, medium value species; LS, low values species.

one-against-one approach for the classification of HS, MS, and LS habitat species classes. Once SVM classification probabilities were calculated, an uncertainty measurement L(x) based on entropy was derived using (12) L(x) = −

K 

pk (x)log2 pk (x)

(12)

k=1

with the only difference, compared with (1), that K = 3, and pk (x) represents the classification probabilities of object x belonging to the class k. A smaller value of L(x) implies that the SVM classification for object x is more reliable. If L(x) value is large, the SVM decision for the object x may be not reliable enough to give a single class label. In this case, the segment is regarded as an unknown segment at level 2. Again, we applied a simple threshold TL2 to separate known and unknown SVM classified objects (i.e., ITCs). Entropy measures the minimum number of bits needed to encode/describe the classes based on the probability of each class. For example, if there are four classes and each class has the same probability, 2 bits (00, 01, 10, 11) is needed. If there is one class with 100% probability of its occurrence, we know the class with full certainty and 0 bits is needed. In our case, three classes with unequal probabilities require less than 2 bits. The range of L is determined as [0 1.59]. In this study, the threshold TL2 = 1.2 was experimentally chosen. If L(x) > TL2 , the segment x was regarded as a level 2 unknown object and given a unique label of UL2 . Otherwise, we assigned x the unique class label that was already determined by SVM. As a result, an object-based classification map at level 2 can be obtained with each segment labeled as one of {“HS”, “MS”, “LS”, and “ UL2 ”}. Accuracy of the SVM classification was assessed by the tenfold cross validation approach. A confusion matrix is reported showing the OA, omission, and commission errors. At this point, each ITC had a class label of either UL1 (i.e., the unknown class at level 1), or “coniferous canopy” at level 1. At level 2, each ITC had a label among the classes: HS, MS, LS, and UL2 (i.e., the unknown class at level 2). The following rules were proposed to make a final decision for each ITC (see Table I). 1) For any object x, if the label was UL1 at level 1 and UL2 at level 2, the current information was insufficient to determine its species type, because the classification at both spatial levels were unknown. Therefore, the object was given a final class label of “unknown.”

Fig. 8. Pixel level classification by ML algorithm.

Fig. 9. Segmentation result of one forest site. The red polygon shows the initial segmentation boundary. The color filled regions indicates the remaining objects after the procedure of removing nonconifers, based on pixel level classification.

2) If the label for the object at level 1 was “coniferous canopy” whereas it was UL2 at level 2, we labeled the object as “unknown conifer.”

1494

IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING, VOL. 8, NO. 4, APRIL 2015

TABLE II SVM C LASSIFICATION ON T RAINING S AMPLES

OA, overall accuracy; HS, high value species; MS, medium value species; LS, low values species.

3) If the label for the object was not UL2 at level 2, we trusted the SVM classification result and labeled the object as the same as that from classification at level 2. In addition, to investigate how important each feature group in term of its contribution to the classification accuracy, we performed three additional SVM classification tests using different types of features. We started with spectral features as the initial classification test and then added one feature group at a time. In each of these tests, we used the grid search again. To assess the final product map, the percentages of correctly classified samples per species group, of unknown conifers, and of unknown species were computed. In addition, 11 circular polygons with 20 m radius (a basic plot size in forest resources inventory) within the study area were manually digitized on imagery, representing areas concentrated by HS, MS, and LS. The composition of each species group class for each polygon was estimated by other trained interpreters, recorded as a percentage value which was rounded to the closest 10%. We then compared our final predicted results with the manual digitized plots. IV. R ESULTS The level 1 classification result is shown in Fig. 8. Based on visual examination and comparison with the original 0.4 m high resolution aerial image, it is shown that most pixels were correctly classified. A few omission errors occurred on an area where a large fraction of small larch and spruce trees were present. The overall classification accuracy was 83% using a total of 500 randomly selected pixels as testing samples. Segmentation result is shown in Fig. 9, where deciduous, grass, and other nonvegetation objects have been masked out. For the SVM classification at level 2, the overall training accuracy of the classification using all 181 training samples was 87% (Table II). The calculated Kappa coefficient was 0.80. Based on tenfold cross validation, the average OA of the ten classification trials was 76.8%. Instead of the training accuracy, the cross validation accuracy was used to represent the general performance of the SVM classifier in our study case. The OA of the cross validation for each classification test for the feature importance investigation is reported in Table III. Table III also shows that the GLCM statistical features and proposed LBPI improved the cross validation accuracy from 72% to 77%. Although this may not be a statistically significant improvement in terms of classification accuracy, the potential of proposed texture features for our classification problem is visible.

TABLE III SVM C LASSIFICATION U SING C OMBINED F EATURE G ROUP

SP, spectral feature; GLCM, gray-level cocurrent matrix feature; LBPI, local binary index feature; OA, overall accuracy.

Fig. 10. Final classification result of one forest site. The red polygon shows the initial segmentation boundary. The color filled regions indicates the class of the remaining objects after the decision fusion procedure.

The obtained final classification map after the procedure of decision fusion is shown in Fig. 10. The second data column in Table IV shows that only 4% of the 181 samples trees were incorrectly treated as nonconiferous objects in the level 1 process, and 3% were identified as unknown objects by the decision fusion process. In addition, 6% of the sample trees were identified as unknown conifers, likely representing those samples hard to be classified. Table V implies that the fraction of predicted species groups (e.g., column HS_p in Table V) is generally in a good agreement with the interpreters’ identification (e.g., HS_i in Table V) except for the result of plot 10. Most of the plots have 0%–10% difference between the predicted and interpreted species groups. It is also observed from Table V that the mean percentage of our predicted HS and LS classes were lower than the interpreted ones by about 5%–10% (48% vs. 38% for HS, 18% vs. 13% for LS). One major reason is that trees that are hard to be classified were likely included into the unknown conifer or unknown object group (7%–8% in average, Table V). It is noted that the visually identified species composition based on image interpretation varies with respect to the interpreter’s skill and imagery quality; therefore, we only intend to use it for a relative comparison with our prediction instead of for the validation of the absolute decision fusion accuracy. If more field data were available, a better validation procedure could be given.

LI et al.: TWO-LEVEL APPROACH FOR SPECIES IDENTIFICATION OF CONIFEROUS TREES

TABLE IV P ERCENTAGE OF S AMPLES P REDICTED BY D ECISION F USION

1495

classification accuracy will be conducted using the greenness index proposed by [10] and adjusting the identified shaded pixels. B. Classification

Unknown C, unknown conifers; HS, high value species; MS, medium value species; LS, low values species. TABLE V C OMPARISON OF P ERCENTAGES OF P REDICTED AND I NTERPRETED S PECIES G ROUPS

HS, high value species, MS, medium value species, LS, low values species; suffix_i indicates results from image interpretation; suffix_p indicates results from prediction by our proposed method.

V. D ISCUSSION A. Segmentation and Imaging Geometry A satisfying ITC segmentation accuracy (73%) has been reported in [21]. As one would expect, the segmentation errors affected the values of the remote sensing features calculated and thus classification result. ITC boundaries may include a few of noncrown pixels due to spectral mixing. This problem can be well solved if high density LiDAR data are available. The imaging geometry may have effects on species identification using the proposed method as well. Near the centers of an image tile (along the cross-track direction), tree crowns are well preserved with little geometric distortion. Tree crowns at the edges of the image tile may lean and shadows may be more apparent. The shaded area within a crown may affect the calculated feature values and cause some misclassification. However, because the solar zenith angle during the image data acquisition was very small, the shaded areas in the image were not significantly large. The proposed method, therefore, works effectively on the given images in this study. Nevertheless, cautions should be taken to apply this method on an image acquired with a large solar zenith angle and thus with a large number of shaded pixels. A possible way to minimize this effect and improve species identification accuracy is to perform true-orthorectification on the given images before using the proposed method. In the future work, an investigation on how the crown delineation affects the feature significances and overall

The main objectives of level 1 classification were to exclude nonvegetation areas and clumps of deciduous trees with a high confidence. The nonvegetation areas were identified very well in the classification because of the differing spectral signature in near infrared region between vegetated and nonvegetated areas. The level 1 classification showed limitations concerning the recognition of isolated deciduous trees featuring a crown sizes smaller than the 4 m pixel size. However, as sparse and isolated deciduous trees do not significantly affect the winter habitat of white-tailed deer from an operational and project scope, we believe that our level 1 classification accuracy satisfied the purpose of modeling deer habitat in wide range ecosystems. As it was noticed that hyperspectral data with moderate pixel resolution (e.g., 4 m in [5]) are becoming available and more and more utilized in the tree species identification research field, such data may performs better than the ADS40 imagery for the first stage classification due to the availability of more spectral bands and enriched spectral information to discriminate between species. Although hyperspectral data alone may be not sufficient for object-based classification at crown-level because of its spatial resolution limit, it may have great potential in improving our proposed framework if combined with submeter high resolution optical imagery. Before hyperspectral data are applied for operational applications at the ITC level, there are two key aspects that require further investigations: 1) data volume and computation costs and 2) effective data fusion methodologies. The decision fusion approach in this study is an open framework that can be easily extended from two-level to multilevel. If additional remote sensing or geographic information such as contextual and elevation layers are available, our fusion rules can be possibly improved using the majority voting approach or more sophisticated fusion theory such as Dempster–Shafer theory. The reason we proposed the identification of “unknown conifer” class in level 2 was under the consideration that users can have their own choice either merging it into the “unknown” class or retaining it as an individual separate class from others. The two thresholds TL1 and TL2 were important for the final classification quality and they had impacts on the uncertainty and the number of remaining objects on the final classification map. At current stage, these two parameters are open to users. Increasing these thresholds may potentially make the classified HS, MS, and LS segments more reliable, whereas the number of those reliable segments may decrease, i.e., there would be more segments labeled as “unknown” category. During the SVM training stage, all training samples were included and labeled as one of the three classes, as a result, there were no unknown objects in the SVM classification validation (i.e., tenfold cross validation). In the final predicted species map, the fraction of unknown objects is 0.11 (i.e., 11%, Table IV) when the threshold TL2 = 1.2 was used. Increasing the threshold

1496

IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING, VOL. 8, NO. 4, APRIL 2015

results in decreasing the faction of unknown objects while the overall classification becomes less reliable. In contrast, decreasing the threshold results in more reliable overall classification and better accuracy but the faction of unknown objects will increase. However, in many research cases, a large number of unknown species is not preferred. A tradeoff between the classification uncertainty and number of reliable segments needs to be determined by users based on their own experiments using their available image data. It is possible that these two parameters could be automatically determined based on statistical analysis of S and L of all segments in the study area. Further research on this will be conducted in our future works. The number of sample trees in the HS class was quite larger than those in median and LS. There may be a potential negative effect that the SVM classification may favor the majority class (i.e., Hemlock). However, this effect seems not very critical in this study as we did not observe obvious differences of the classification accuracy between the majority and minority classes. In addition, this potential effect could be reduced by applying over- or down-sampling techniques such as the synthetic minority over-sampling technique (SMOET). However, whether or not the classification accuracy will be decreased requires more significant tests using these different sampling techniques. Keeping in mind that different remote sensing data and methods, species, and site conditions from existing studies limit direct comparability, we achieved similar classification accuracy as some of other studies using optical imagery (e.g., [11], [13], and [31]). In Ontario, wildlife habitat modeling mostly relies on land cover classification map and forest resources inventory, in which detailed species information of individual trees are unavailable. Our study provides an initialization toward the ecology modeling using high spatial resolution aerial imagery and it is likely that this attempt will improve the modeling performance. Our method can be directly transferred and used for large forest areas and operational applications; however, more evaluation procedures should be conducted before transferring. VI. C ONCLUSION Forest resources inventory and wildlife habitat conservation require detailed and accurate information on species or species groups. We proposed an effective and flexible framework which employs information from both low and high spatial resolution imagery to identify coniferous species groups that will contribute to wildlife habitat modeling and management. The decision of classification at both pixel and object levels was combined based on measures of entropy and simple logic rules. The structure-based LBPI feature was designed and proved to be effective for improving species classification accuracy. The success of this research delivers an alternative way to obtain species information of individual trees when airborne LiDAR data are not available. ACKNOWLEDGMENT The authors would like to thank B. Naylor and L. Landriault for their contributions in identifying the characteristics of the winter habitat white-tailed deer, and B. Naylor for providing the training and testing samples.

R EFERENCES [1] T. Allouis, S. Durrieu, C. Vega, and P. Couteron, “Stem volume and above-ground biomass estimation of individual pine trees from LiDAR data: Contribution of full-waveform signals,” IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., vol. 6, no. 2, pp. 924–934, Apr. 2013. [2] B. Hu, J. R. Miller, P. Zarco-Tejada, J. Freemantle, and H. Zwick, “Boreal forest mapping at the BOREAS study area using seasonal optical indices sensitive to plant pigment content,” Can. J. Remote Sens., vol. 34, pp. s158–s171, 2008. [3] M. Dalponte et al., “Tree species classification in the Southern Alps based on the fusion of very high geometrical resolution multispectral/ hyperspectral imaged and LiDAR data,” Remote Sens. Environ., vol. 123, pp. 258–270, 2012. [4] M. A. Cho et al., “Mapping tree species composition in South African savannas using an integrated airborne spectral and LiDAR system,” Remote Sens. Environ., vol. 125, pp. 214–226, 2012. [5] A. Ghosh et al., “A framework for mapping tree species combining hyperspectral and LiDAR data: Role of selected classifiers and sensor across three spatial scales,” Int. J. Appl. Earth Observ. Geoinf., vol. 26, pp. 49–63, 2014. [6] R. Dinuls, G. Erins, A. Lorencs, I. Mednieks, and J. Sinica-Sinavskis, “Tree species identification in mixed baltic forest using LiDAR and multispectral data,” IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., vol. 5, no. 2, pp. 594–603, Apr. 2012. [7] R. Palenichka, F. Doyon, A. Lakhssassi, and M. Zaremba, “Hierarchical multi-scale segmentation of LiDAR images in forest areas,” in Proc. IEEE Int. Geosci. Remote Sens. Symp. (IGARSS), Jul. 22–27, 2012, pp. 5462– 5465. [8] U. Benz, P. Hofmann, G. Willhauck, I. Lingenfelder, and M. Heynen, “Multi-resolution, object-oriented fuzzy analysis of remote sensing data for GIS-ready information,” ISPRS J. Photogramm. Remote Sens., vol. 58, pp. 239–258, 2004. [9] J. Li, B. Hu, M. Woods, B. Naylor, and L. Landriault, “Integrating spectral and textural features derived from high resolution aerial imagery for species classification of individual trees,” presented at the 34th Can. Symp. Remote Sens., Victoria, BC, Canada, Aug. 27–29, 2013. [10] K. Zhang and B. Hu, “Individual urban tree species classification using very high spatial resolution airborne multi-spectral imagery using longitudinal profiles,” Remote Sens., vol. 4, pp. 1741–1757, 2012. [11] D. Leckie et al., “Automated tree recognition in old growth conifer stands with high resolution digital imagery,” Remote Sens. Environ., vol. 94, pp. 311–326, 2005. [12] T. Brandtberg, “Classifying individual tree species under leaf-off and leaf-on conditions using airborne LiDAR,” ISPRS J. Photogramm. Remote Sens., vol. 61, no. 5, pp. 325–340, 2007. [13] L. T. Waser, C. Ginzler, M. Kuechler, E. Baltsavias, and I. Hurni, “Semiautomatic classification of tree species in different forest ecosystems by spectral and geometric variables derived from airborne digital sensor (ADS40) and RC30 data,” Remote Sens. Environ., vol. 115, no. 1, pp. 76–85, 2011. [14] M. Immitzer, C. Atzberger, and T. Koukal, “Tree species classification with random forest using very high spatial resolution 8-band WorldView2 satellite data,” Remote Sens., vol. 4, pp. 2661–2693, 2012. [15] T. A. Warner, J. Yeol Lee, and J. B. McGraw, “Delineation and identification of individual trees in the eastern deciduous forest, in Proc. Automated Interpretation of High Spatial Resolution Digital Imagery for Forestry, Canadian Forest Service, Pacific Forestry Center, Victoria, BC, Canada, 1998, pp. 81–91. [16] S. E. Franklin, Remote Sensing for Sustainable Forest Management, 1st ed. Boca Raton, FL, USA: CRC Press, 2001. [17] C. A. Coburn and A. C. B. Roberts, “A multiscale texture analysis procedure for improve forest stand classification,” Int. J. Remote Sens., vol. 25, no. 20, pp. 4287–4308, 2004. [18] M. S. Kulikova, M. Mani, A. Srivastava, X. Descombes, and J. Zerubia, “Tree species classification using radiometry, texture and shape based features,” in Proc. Eur. Signal Process. Conf. (EUSIPCO), 2007, pp. 1–5. [19] Z. Li, R. Hayward, Y. Liu, and R. Walker, “Spectral–texture feature extraction using statistical moments with application to object-based vegetation species classification,” Int. J. Image Data Fusion, vol. 2, no. 4, pp. 347–361, 2011. [20] J. Li, B. Hu, and T. Noland, “Classification of tree species based on internal structure features derived from high density LiDAR data,” Agric. Forest Meteorol., vol. 171–172, pp. 104–114, 2013. [21] L. Jing, B. Hu, T. Noland, and J. Li, “An individual tree crown delineation method based on multi-scale segmentation of imagery,” ISPRS J. Photogramm. Remote Sens., vol. 70, pp. 88–98, 2012.

LI et al.: TWO-LEVEL APPROACH FOR SPECIES IDENTIFICATION OF CONIFEROUS TREES

[22] R. M. Haralick, K. Shanmuga, and I. Dinstein, “Textural features for image classification,” IEEE Trans. Syst. Man Cybern., vol. 3, no. 6, pp. 610–621, Nov. 1973. [23] Z. H. Guo, L. Zhang, and D. Zhang, “Rotation invariant texture classification using LBP variance (LBPV) with global matching,” Pattern Recognit., vol. 43, no. 3, pp. 706–719, 2010. [24] T. Ojala, M. Pietikainen, and T. Maenpaa, “Multiresolution gray-scale and rotation invariant texture classification with local binary patterns,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 24, no. 7, pp. 971–987, Jul. 2002. [25] H. Mills, “Analysis of the transferability of support vector machines for vegetation classification,” in Proc. Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., 2008, vol. 37, pp. 557–563. [26] H. Lin, C. Lin, and R. Weng. (2001). A Note on Platt’s Probabilistic Outputs for Support Vector Machines [Online]. Available: http://www.csie.ntu.edu.tw/~htlin/paper/doc/plattprob.pdf, accessed on Oct. 23, 2014. [27] J. Platt, “Probabilistic outputs for support vector machines and comparison to regularized likelihood methods,” in Advances in Large Margin Classifiers, A. Smola, P. Bartlett, B. Scholkopf and D. Schuurmans, Eds. Cambridge, MA, USA: MIT Press, 2000. [28] A. Karatzoglou, D. Meyer, and K. Hornik, “Support vector machine in R,” J. Stat. Softw., vol. 15, no. 9, 2006, pp. 1–28. [29] T. Wu, C. Lin, and R. Weng, “Probability estimates for multi-class classification by pairwise coupling,” Adv. Neural Inf. Process., vol. 16, pp. 975–1005, 2003. [30] F. Melgani and L. Bruzzone, “Classification of hyperspectral remote sensing images with support vector machines,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 8, pp. 1778–2004, Aug. 2004. [31] B. Mora, M. Wulder, and J. White, “Identifying leading species using tree crown metrics derived from very high spatial resolution imagery in a boreal forest environment,” Can. J. Remote Sens., vol. 36, no. 4, pp. 332– 344, 2010.

Jili Li received the B. Eng. and M. Eng. degrees in remote sensing and photogrammetry from Wuhan University, Wuhan, China, in 2005 and 2008 respectively, and the Ph.D. degree in earth and space science from York University, Toronto, ON, Canada, in 2013. He worked as Research Assistant and Teaching Assistant at York University from 2008 to 2013 and Postdoc Fellow at 2014. He is currently an NSERC Industrial Postdoc Fellow with KBM Resources Group, Thunder Bay, ON, Canada. His research interests include LiDAR and multispectral remote sensing for precision forest resources inventory, image processing, and machine learning.

1497

Baoxin Hu (M’03) received the Bachelor and Master degrees in electrical engineering from Tianjin University, Tianjin, China, in 1987 and 1990 respectively, and the Ph.D. degree in remote sensing and GIS from Boston University, Boston, MA, USA, in 1998. She is currently an Associate Professor with the Geomatics Engineering, York University, Toronto, ON, Canada. Her research interests include individual tree species delineation and classification and the retrieval of biophysical and biochemical parameters of vegetation canopies from remotely sensed data, and 3-D scene reconstruction from LiDAR and optical imagery.

Murray Woods received the degree in forest technologist from Sir Sandford Fleming College, Peterborough, ON, Canada, in 1983. He is a Senior Analyst with the Ontario Ministry of Natural Resources, ON, Canada. He started his career with the Ontario Forest Research Institute participating in growth and yield modeling for Ontario’s Southern, Central, and Boreal forest regions. In 2006, he assumed his current role where he is working with a team exploring the utility of LiDAR technology and multiband digital image analysis in the development of enhanced forest inventories at the stand and at the individual tree level.