MAPPING FOREST STAND COMPLEXITY FOR WOODLAND CARIBOU HABITAT ASSESSMENT USING MULTISPECTRAL AIRBORNE IMAGERY

The decline of the woodland caribou population is a result of their habitat loss. To conserve the habitat of the woodland caribou and protect it from extinction, it is critical to accurately characterize and monitor its habitat. Conventionally, products derived from low to medium spatial resolution remote sensing data, such as land cover classification and vegetation indices are used for wildlife habitat assessment. These products fail to provide information on the structure complexities of forest canopies which reflect important characteristics of caribou’s habitats. Recent studies have employed the LiDAR system (Light Detection And Ranging) to directly retrieve the three dimensional forest attributes. Although promising results have been achieved, the acquisition cost of LiDAR data is very high. In this study, utilizing the very high spatial resolution imagery in characterizing the structural development the of forest canopies was exploited. A stand based image texture analysis was performed to predict forest succession stages. The results were demonstrated to be consistent with those derived from LiDAR data. * Corresponding author.


INTRODUCTION
Woodland caribou (Ranifer Tarandus Caribou) is an important wildlife species ecologically and culturally in Canada.This amazing creature is now at risk and listed as a threatened species provincially and nationally (Ontario Woodland Caribou Recovery Team, 2008).The decline in woodland caribou populations is primarily the result of habitat loss and forest fragmentation which is caused by human interference including spread of agriculture, oil and gas exploration, and mining (McLoughlin, 2003).To protect them from the threat of extinction, it is critical to conserve its habitat conservation.
To map the wildlife habitat, the basic idea is to establish a link between the organism's characteristics and behaviors to the physical habitat (Guisan and Zimmermann 2000).Conventionally, the observation was conducted by field work which has the promising accuracy and generally been used on the local scale.For the past decades, remote sensing techniques played a key role in wildlife habitats mapping, when large spatial information becomes more favourable for broader area analysis while conventional field work is time and labour intensity.(Kerr and Dostoevsky, 2003).The forest inventory map and land use map which generated from lower resolution satellite imagery has been employed as the major data source to map the woodland caribou habitat (Hansen, 2001).In addition to land cover/use maps, other remote sensing products, such as NDVI (Normalized Difference Vegetation Index) images calculated from satellite data with medium to low spatial resolution have been used in other wildlife habitat mapping.However, the map was not usually up to date, and classes from products usually generated for other or generic applications are not best suited for woodland caribou habitat mapping.These products are not sufficient in characterizing caribou or other wildlife habitat at local (fine spatial resolution) scales and in characterizing sub-canopy vegetation structure relevant to wildlife (Brown, 2006) Forest vertical structure influences animal-habitat associations and biodiversity significantly (Brokaw and Lent 1999).Recently, studies have focused on the exploitation of aerial digital imagery and high spatial resolution satellite imagery in characterization of forest structure.Song and Woodcock (2002) characterized the sizes of tree crown using the sill of a variogram.Pasher and King (2010) used multivariate texture analysis to produce a continuous map of the structural complexity, which was based on the empirical relationship between the texture variables and field measurements.Kayitakire et. al (2006) used 1-m resolution IKONOS-2 imagery to estimate the forest variables from the texture analysis of Grey Level Co-occurrence Matrix (GLCM).Their results showed that the structure characteristics including top height, circumference, stand density and age variables are correlated to texture variables.Ouma et al (2008) employed grey-level cooccurrence matrix (GLCM) and wavelet transform (WT) texture analysis for the differentiation of forest and non-forest vegetation type using QuickBird imagery.Their results showed that with the combination of GLCM-mean texture (microtextures) WT-derived texture (macrotextures), the differentiation and classification of the overall vegetation types can be improved.
At longer temporal scale, forest succession is more primary than the structure, as the structure only capture the forest at a certain time period while the forest succession describe the forest development dynamics (Shugart, 2000).In addition, accurate classifications of forest succession stage at large spatial extents are critical to characterize wildlife habitat (Helle and Monkkonen, 1990).Current Studies have shown that a combined set of LiDAR derived height metrics which characterizing the three-dimensional structure of forest canopies is able to map forest successional stage with a high accuracy (Falkowski, 2009;van Ewijk, 2011).
Although LiDAR data have the ability to reveal the succession stage of the forest stands, the data acquisition and processing take greater cost (Jensen 2005).High resolution remote sensing imagery comes at a greater availability, but no study has explored the relationship in between the image texture and forest succession stages.The objective of this study is to classify forest succession stage through very high resolution multispectral aerial imagery.Specifically, image texture variables which were derived from Gray Level Co-occurrence Matrix and shadow fraction were employed to classify the forest succession stages.In addition, to determine the reliability of the result derived from the image texture analysis, cross validation with LiDAR data was performed.Existed LiDAR metrics was adopted to predict the forest succession stages, and then compared with the classification result generated from imagery for accuracy assessment.

STUDY AREA AND DATA SET
The study area is located near the town of Hearst, Central Ontario, Canada (49.7°N, 83.7°W).Hearst forest falls within the boreal mixed wood region and covers approximately 1.23 million ha; 1.00 million ha of which is productive forest.The study site contains of black spruce (picea mariana), white spruce (picea glauca), balsam fir (abies balsamea), trembling aspen (populus tremuloides) and black poplar (populus balsamifera).The optical imagery which was acquired using a Leica ADS-40 in the summer of 2007 has four spectral bands, blue, green, red and near infrared with a spatial resolution of 0.4 m by 0.4 m.LiDAR data were collected using an Optech ALS50 sensor mounted in a Cessna 310 aircraft in summer 2007 during leaf-on conditions.The LiDAR data were discrete pulse, with an average 1.677 pulses per square meter.Location of the study area is shown in Figure 1 below with a snapshot of the multispectral imagery.

METHODOLOGY
To differentiate the stages in forest development, a suite of variables were extracted from the imagery representing the variance within-crown, within-shadow, and canopy level spectral and spatial properties.According to the literature, two approaches have been popular to use, which are semivariograms analysis and Gray Level Co-occurrence Measure (GLCM).Semivariogram contain measurements that relates to forest structure (Lé vesque and King, 2003;Treitz and Howarth, 2000).The range of the semivariogram is an indicator of the distance of spatial dependence, often related to the size or scale of dominant objects in the imagery (Curran, 1988;Johansen et al., 2007), and the sill was proved to be an indicator of the total structurally dependent variance in the data that was expected to respond in a manner similar to the texture measures (Johansen et al., 2007).However, semivariogram calculation is complex and computation extensive, which is not suitable for large area mapping.Therefore, the semivariograms for imagery over the whole study area for forest succession stages were not calculated.GLCM approach is adopted for this study.The GLCM variables were calculated first, and used for classifying different forest developing stages.
Shadow fraction, which indicates the tree canopy closure and tree densities, is also an important variable to measure the forest structure.Therefore, shadow fraction was calculated and introduced as an additional variable to the forest development stages classification to test the accuracy improvements.Details of variables extraction and analysis are given in the following paragraphs.
In this study, four general stage of the forest development was adopted.The stand initiation stage indicates where new vegetation becomes established and fully occupies the disturbed site.The young multi-storey stage follows and is mainly characterized by competition among the dominant trees.After the time goes by, young trees start to regenerate and the stand enters the understory re-initiation stage.In the last stage, state as the old growth stage, autogenic and gap-replacing processes create patches large enough for stand initiation to begin (Oliver and Larson, 1996).
From the entire study area, 20 test areas representing four different succession stages, with the size of 501 x 501 pixels (200m x200m) were subset based on visual interpretation.Example of each forest succession stages are listed in Figure 2. In this study, a size of 400 square meters is used to represent a stand.To capture the spectral and spatial variation within a stand, all 20 test areas are arbitrarily cut into numbers of blocks with a size of 50 pixels by 50 pixels and processed for the image texture analysis.

GLCM Calculations
For very high resolution aerial imagery, the gaps and shadings in between tree crowns are visible and shown as the spatial variability in image brightness.As the forest developed from young to mature, the amount of the gap and shadings changed as well.Therefore, texture measurements were performed through statistical calculated GLCM.As reviewed from the literature, homogeneity was more effective than the first-order variance in discriminating age classes of Douglas-fir forest stands from IKONOS panchromatic data (Franklin et al., 2001), others also found that the contrast , ASM, entropy and homogeneity measures should be used for forest attributes.(Pesaresi 2000, Cosmopoulos and King, 2004, Kayitakire et al., 2006).
Three parameters controlling the statistics calculated from GLCM were considered in this study, moving window size, displacement and the direction angle.Various values of these parameters have been considered to cover the possible range, while maintaining a manageable number of texture cases.The four main directions (0°, 45°, 90° and 135°) and window sizes 5pixels ×5 pixels, which was intend to capture the variation within crowns.The displacement values were set to 1.The eight texture measures were computed and then further averaged over a plot of 400 square meters, which is considered as the stand size in this study.
Where is the th entry of the normalized GLCM matrix; , , and are the mean the standard deviation respectively.
To investigate differentiating all four forest succession stage from the image texture variables, random forest classification was employed.The random forest classifier consists of a combination of tree classifiers where each classifier is generated using a random vector sampled independently from the input vector, and each tree casts a unit vote for the most popular class to classify an input vector (Breiman 2001).Recent applications have shown that random forest is good for classifying hyper spectral remote sensing data, and multisource geographic data (Ham et., al., 2005, Gislason et al., 2006, Lawrence et. al., 2006)

Shadow fractions
Different from the low resolution satellite imagery, shadows can be observed in high resolution aerial imagery.The amount of the shadow varied when the tree height, presence of understory and size of the overstory are different.The shadow fraction can be used as one descriptor to measure the structure difference of the forest.In our study, shadow fraction was calculated as the ratio of shaded area to total area within the window size 50 x 50 pixels (20m x 20m).It was derived using the NIR band because it contained the highest contrast between shaded and nonshaded forest areas among all four spectral bands.
Based on the distribution of the histogram, a threshold value was defined to separate the pixels into shades and non-shades groups.The ration between the number of shaded pixel and the total number of pixels within the window is the shadow fraction.
Averaged results are shown in Table 2.
As shown in the table 2, each stage of the forest shows different percentage of shadow presences, while the difference in between forest stages is not significant to separate from one another.To test the importance of this texture variable, the shadow fraction was introduced as an additional variable into the classification.Together with the variables derived from the GLCM, the random forest classification accuracy was improved to 92.18%.Also, among all the variables, the shadow fraction was ranked the third important variable in random forest classification.Overall, the shadow fraction does increase the classification.Specifically, adding the shadow fraction improved differentiate the old growth stage and understory reinitiation, and the stand initiation and young multi-storey.However, it does not improve classify the young multistory and understory re-initation.The confusion matrix of the classification result is shown in

RESULT VALIDATION
To validate the result from the image texture measurements, a cross validation was applied against the LiDAR data.Since LiDAR data measure the three-dimensional forest structures, such ability makes it possible to classify the forest succession stages.Figure 3 shows the extracted LiDAR points and height distribution.

CONCLUSION
This study shows that using GLCM texture measurements alone can differentiate stages of forest stands development with achieved classification accuracy over 89%.The reason account for that is, as the forest stands develop, the canopy surface will change from smooth to rough, which makes the variation in between canopies measureable from the textural analysis.Specific correlation is sensitive to capture the dependency, thus able to distinguish the homogenous area and the areas with greater variation.In addition, shadow fraction can measure the amount of tree and also the density of the tree crowns.It can be used roughly describe the forest development, but can not distinguish all four stages clearly.However, employed as an additional variable to the variables derived from GLCM, classification result accuracy can be improved to 92%.As the forest developed from young multistory to understory reinitation, the understories are still too small to reflect the sunlight, thus can be taken into account as shadows.Therefore, adding the shadow fraction cannot improve the classification in between these two stages.
In addition, as the LiDAR data has the ability of viewing the forest stands vertically.It was used to statistically describe the development of the forest.And to further test the accuracy of the from the texture analysis, cross validation was performed from the classification result and the LiDAR derived indices.As the class from young to mature corresponding to the predicted Lorey's height from low to height, the result derived from texture analysis agreed with the LiDAR derived indices.Therefore, the ability of using texture measurements of the high spatial resolution imagery to derive forest complexity information was demonstrated.
Future work will include analyze how GLCM parameters, interpixel displacement and the window size influence the estimates of forest structures.At different stages of the forest, the canopy openness, size of the crow/stands will be different.At the early stage of the forest, the trees are more isolated, while in the mature stage, the trees grow into each other, and the overlaps of the tree crown make multiple crowns clustered as one single object.Thus, one window size can not fit all forest stages.Although in this preliminary analysis, a standard window size provides promising result, further analysis will include segmentation process to extract canopy objects.The texture measurements will be processed within each object rather than predefined windows.

Figure 1 .
Figure 1.Research area of Hearst Forest, Ontario, Canada

Figure 2 :
Figure 2: Description of the four forest stand development stages.From top to down and left to right: (a) stand reinitiation, (b)stem exclusion, (c) understory re-initiation, and (d) old growth stage.The false colour composite of the optical imagery covering the study area with the near-infrared band displayed as red, red as green and green as blue

Figure 3 :
Figure 3: Discrete return LiDAR points with height distribution for four forest succession stages.From top to down and left to right: (a) stand initiation, (b)young multistory, (c) understory reinitiation, and (d) old growth stage.

.
Eight texture variables, Mean, Variance, Contrast, Entropy, ASM, Homogeneity, Dissimilarity, and Correlation, at four direction angles are used together as the 32 parameters in random forest.After the classification, 1433 out of 1600 samples are correctly classified, which is in a result of 89.56 % accuracy.The confusion matrix of the classification result is shown in Table1.

Table 1 .
Confusion Matrix from random forest classification

Table 3 :
van Ewijk et al (2011)mpirical model should be built to characterize the forest succession stage.Due to the lack to filed plot, no sufficient information was available to establish a statistical model of using LiDAR data representing forest stage.The model (Equation (1)) reported invan Ewijk et al (2011)to predict Lorey's height for forest succession stage was adopted.Based on their findings, the Lorey's height which is basically the stand tree height weighted by basal area was found to be a good index to separate the young multistory stage from an understory reinitiation stage.To validate the forest succession stage derived from the texture measurements, 400 samples are employed for the classification from random forest using the previous defined training data.On the other hand, the LiDAR points from the overlapped area from the samples are used to calculate the predicted Lorey's height.Averaged result is shown in Table4below.The results demonstrate that the classification results and LiDAR derived predicted Lorey's height agree with each other.

Table 4 .
Cross validation from the texture measurement classification result with the LiDAR derived predicted Lorey's height.