Mapping Savanna Tree Species at Ecosystem Scales Using Support Vector Machine Classification and BRDF Correction on Airborne Hyperspectral and LiDAR Data

Abstract

Mapping the spatial distribution of plant species in savannas provides insight into the roles of competition, fire, herbivory, soils and climate in maintaining the biodiversity of these ecosystems. This study focuses on the challenges facing large-scale species mapping using a fusion of Light Detection and Ranging (LiDAR) and hyperspectral imagery. Here we build upon previous work on airborne species detection by using a two-stage support vector machine (SVM) classifier to first predict species from hyperspectral data at the pixel scale. Tree crowns are segmented from the lidar imagery such that crown-level information, such as maximum tree height, can then be combined with the pixel-level species probabilities to predict the species of each tree. An overall prediction accuracy of 76% was achieved for 15 species. We also show that bidirectional reflectance distribution (BRDF) effects caused by anisotropic scattering properties of savanna vegetation can result in flight line artifacts evident in species probability maps, yet these can be largely mitigated by applying a semi-empirical BRDF model to the hyperspectral data. We find that confronting these three challenges—reflectance anisotropy, integration of pixel- and crown-level data, and crown delineation over large areas—enables species mapping at ecosystem scales for monitoring biodiversity and ecosystem function.

1. Introduction

Savannas harbor spatially complex assemblages of vegetation that are mediated by an array of biotic and abiotic factors including plant competition, fire, herbivory, soils and climate [1–4]. Mapping the distribution of species abundances across spatial scales relevant to these processes is requisite to understanding their role in shaping and maintaining savanna biodiversity. Remote sensing continues to increase in importance in providing the spatial information needed to quantify and monitor spatial and temporal variation in savanna vegetation, but few approaches have yielded detailed data on species occurrence and distribution. Nevertheless, high spatial resolution airborne imaging spectrometers and Light Detection and Ranging (LiDAR) sensors are promising in this regard because they allow detailed spectral and structural analysis of vegetation at large scales [5–9].

Airborne imaging spectroscopy has recently been used to discriminate different tree species, which can serve as the starting point for mapping spatially explicit patterns of biodiversity in the context of environmental controls. Cho et al.[6] used the spectral angle mapper (SAM) classifier to minimize intraspecific spectral variability to discriminate among South African savanna tree species, while Naidoo et al.[10] integrated hyperspectral and LiDAR measurements using the Random Forests approach to improve savanna species identification. Elsewhere, the support vector machine (SVM) classifier has been applied to hyperspectral and LiDAR data to predict tree species in deciduous [7] and tropical forests [11]. The SVM classification technique has been increasingly applied to classification of airborne hyperspectral imagery [12–14], where its higher accuracy compared to traditional techniques (e.g., maximum likelihood, neural networks, decision tree classifiers) stems from a lower sensitivity to high dimensionality (i.e., the Hughes effect) [15]. These and other studies clearly indicate the value of SVM classification and of hyperspectral and LiDAR measurements, yet three issues continue to impact the accuracy of airborne species mapping at large ecosystem scales (e.g., >10,000 ha): anisotropy of land surface reflectance, delineating tree crowns over large areas, and integration of pixel-level with crown-level data.

The reflectance properties of tree canopies vary with view angle and sun position (anisotropy), and the mathematical expression of these properties is known as the bidirectional reflectance distribution function (BRDF) [16–20]. BRDF effects can either be treated as noise or as additional information, depending on the application [21–23]. In the context of species or land cover classification from spectral imagery, BRDF effects are often considered noise due to the variation in reflectance across an image (i.e., cross-track brightness gradients), even if the cover type is constant. This is a particularly important issue with airborne data when multiple flight lines are mosaicked into a single scene, which is often necessary for mapping at larger spatial scales. BRDF models aim to correct for this effect by normalizing each reflectance sample (pixel) to a common view and solar geometry by modeling the light scattering and shadowing occurring within each sample. For example, MODIS multispectral and albedo products are first corrected using the Ross-Li model, a geometric-optical, kernel-based approach [24,25]. Yet the effects of BRDF on the accuracy of species predictions derived from hyperspectral data are unknown, since few studies have applied such anisotropic BRDF corrections to airborne hyperspectral imagery [26–28] due to the additional viewing and solar angle data required for each reflectance pixel and the added complexity of calibrating the kernel coefficients.

Another challenge to mapping species at ecosystem scales is the segmentation and extraction of thousands or potentially millions of polygons, or tree crown objects, from imagery. Previous species classification studies report higher accuracies at the crown-level compared to the pixel-level [11,29], yet training/validation crowns are typically manually delineated, making mapping over large geographic areas impractical. Automated crown segmentation has the potential to facilitate species mapping on a crown-by-crown basis at ecosystem scales over tens of thousands of hectares. However, tree segmentation using spectral data can be problematic in savannas because small trees can have similar LAI to that of shrubs or dense herbaceous plants [30,31]. This distinction also assumes a clear definition of shrub vs. tree; although this is primarily a matter of ecological definition, of which there are several (here we define a tree as greater than 1.5 m in height and greater than 5 cm in stem diameter). LiDAR-based segmentation instead uses vegetation height to detect crown edges and is thus far less susceptible to LAI-related issues. Moreover, LiDAR segmentation allows for extraction of crown-level structural information, such as tree height, which is beneficial for separating species in different size classes. The typical constraint to LiDAR segmentation is that LiDAR data is rarely simultaneously collected with hyperspectral data due to cost and complexity.

Yet integration of crown-level structural characteristics with pixel-level species detection presents a third challenge for traditional classification methods. Previous studies integrating hyperspectral and LiDAR data have used pixel height along with spectral data for species classification [10,32]. Yet the pixel size of airborne data is typically smaller than most tree crowns, leading to pixel heights varying greatly within a tree canopy. Maximum height of the crown (MH) can instead be used to mitigate pixel-level variance in height and to discern shrubs from trees. Crown area (A), volume, and related allometric quantities may also be useful for discriminating among species, yet they have no pixel-level analog.

Here we build upon previous savanna species classification studies in the following ways: (1) a stacked SVM model is used to integrate pixel-based species classification with crown-level LiDAR structural data; (2) an analysis of BRDF effects on species predictions results in a correction that enables ecosystem-scale species mapping across many airborne mapping flight lines collected over multiple days; and (3) an automated crown segmentation of airborne LiDAR data is used to delineate small, low-LAI trees as well as trees orders of magnitude larger. Our overarching goal is to integrate these three advances to achieve higher species mapping accuracies at ecosystem and regional scales.

2. Methods

2.1. Overview

Our proposed processing flow for predicting and mapping species is shown in Figure 1. The methods section provides details for each step, and also for model calibration and several tests not illustrated in the flowchart, namely a comparison of SVM models and spatial effects of BRDF correction on species prediction.

2.2. Study Site

Kruger National Park (KNP) (24°S, 31.5°E) is located in eastern South Africa, spanning an area 360 km north-south and 70 km east-west (Figure 2). KNP is roughly equally split into granite substrates in the west and basalt substrates in the east, with granites weathered to sandy, nutrient-poor soils and basalts weathered to clay-rich, primarily smectitic soils [33]. The climate of KNP is mostly semi-arid with mean annual temperature and precipitation of 22 °C and 550 mm·year⁻¹, respectively, and an average potential evaporation of 7 mm·day⁻¹[2]. Precipitation ranges from ∼350 to 500 mm·year⁻¹ in the north to 500–700 mm·year⁻¹ in the south, with a similar east-west gradient in southern KNP. Five study areas were flown to cover KNP’s major vegetation types, geology, topography, and climatic gradients. Field inventory data of individual trees were then collected within each area to inform and validate airborne species prediction (see Field Data Collection below).

2.3. Airborne Data Collection

The Carnegie Airborne Observatory (CAO) Alpha system [2] collected discrete-return LiDAR data and hyperspectral imagery in April–May 2008. For the LiDAR a pulse repetition frequency of 50 kHz was used to generate three-dimensional maps of tree canopy structure at 1.12 m laser spot spacing (aircraft altitude = 2,000 m). Flights were planned with 100% repeat coverage, and therefore LiDAR point density averaged two points per spot. LiDAR spatial errors were less than 0.20 m vertically and 0.36 m horizontally [34,35]. A woody canopy height map was computed by subtracting a ground DEM (classified from LiDAR last return elevation) from a canopy surface DEM (first return elevation). A physical model was used to calculate both DEMs using the Terrascan/Terramatch (Terrasolid Ltd., Jyväskylä, Finland) software package. Ancillary information necessary for BRDF correction (i.e., view zenith angle, solar zenith angle, and relative azimuth angle) were calculated for each pixel using the LiDAR data and an onboard GPS-IMU (Global Positioning System-Inertial Measurement Unit).

The hyperspectral data were collected using a modified version of the Compact Airborne Spectrographic Imager (CASI-1500), with custom anti-reflective lens coatings to decrease stray light and to boost SNR, a high-throughput read-out to lower integration times, and a VNIR cooling system to increase sensor stability [34]. The sensor consisted of 72 bands from 385 to 1,054 nm at 1.12 m spatial resolution (co-aligned with the LiDAR). Here we define the visible portion of the spectrum (VIS) as 450–650 nm and the near-infrared portion (NIR) as 750–1,050 nm. Surface reflectance was derived from the radiance data using the atmospheric correction software, ACORN 5 BatchLi (Imspec LLC, Palmdale, CA, USA). ACORN used a MODTRAN look-up table to correct for Rayleigh scattering and aerosol optical thickness.

2.4. Field Data Collection

Field inventory data were collected in 2009 within the five overflight areas to inform and validate airborne species identification. Individual tree crowns (n = 729,

Table 1. Summary of field data collected to calibrate and validate the airborne species prediction models. From left to right: woody plant species, family, number of field-identified crowns, total sample size (in pixels) of airborne data (hyperspectral and LiDAR) available per species, and min-max ranges of field-measured basal diameter (cm), crown diameter (m), and height (m).

Species	Family	Tree Crowns	Pixels	Basal D (cm)	Crown D (m)	H (m)
Acacia nigrescens	Fabaceae	45	806	8–72	2–17	3–34
Acacia tortilis	Fabaceae	38	564	5–57	3–17	2–13
Combretum apiculatum	Combretaceae	60	365	2–49	2–11	2–9
Combretum collinum	Combretaceae	29	174	2–24	2–8	1–6
Combretum hereoense	Combretaceae	36	346	2–32	1–7	1–8
Combretum imberbe	Combretaceae	65	3,693	7–135	3–28	5–22
Colophospermum mopane	Fabaceae	44	542	5–96	1–13	2–18
Croton megalobotrys	Euphorbiaceae	12	114	7–90	2–9	2–6
Diospyros mespiliformis	Ebenaceae	31	1,425	16–182	4–28	4–26
Euclea divinorum	Ebenaceae	50	587	1–37	2–9	1–6
Philenoptera violacea	Fabaceae	44	1,397	6–109	2–18	0–19
Spirostachys africana	Euphorbiaceae	26	618	5–53	4–15	2–12
Salvadora australis	Salvadoraceae	26	298	3–45	2–14	1–7
Sclerocarya birrea	Anacardiaceae	73	1,469	4–87	4–20	6–15
Terminalia sericea	Combretaceae	48	408	1–43	2–13	2–11
Other		106	1,306
Total		729	13,998

) were identified for species, located with a hand-held GPS unit (GS50 Leica Geosystems Inc., Norcross, GA, USA), and measured for basal diameter (cm), crown diameter (m), and height (m). The target species shown in Table 1 were collected across as many sites as possible, with most species present at two or more sites. An additional 124 circular field plots of 30 m diameter were used to measure the abundance of each species, but these plots were not used in calibration or validation of the species classifiers since trees were not individually located within the plots.

The species and GPS locations were the only field data used to train the airborne species classifiers, although height and the other field structural measurements demonstrate the maximum size class reached by each species in the calibration/validation dataset (Table 1). Spectral data were extracted from the field-identified crown polygons and filtered to select pixels with NDVI ≥ 0.5 and mean NIR reflectance ≥ 20%. The NDVI filter was used to exclude soil and grass, and the NIR filter excluded heavily shaded samples. Only crowns with three pixels passing both filters were used as model calibration data, which are shown in Table 1. The resulting set of crowns was then grouped into species classes with a minimum of 25 crowns each. Crowns not belonging to one of these species were grouped together into a class labeled “other”. The two exceptions were Acacia burkei and Acacia nigrescens Lannea schweinfurthii, which were grouped with and Sclerocarya birrea, respectively, instead of “other”, due to similar foliar characteristics and canopy structure.

2.5. Crown Segmentation

Crown segmentation was performed using the eCognition (Definiens Developer 8.7) software package, with the following customized “region growing” algorithm and a LiDAR top-of-canopy, vegetation height map as input. In the first segmentation step a pixel was classified as a “crown seed” if both of the following conditions were met: (1) its height > 0.5 m; (2) the surrounding pixels were of similar height. The second criterion was satisfied if the height coefficient of variation (CV) of neighboring pixels was below a given threshold (CV < 0.3). CV was computed for each pixel as the standard deviation of the height of the eight neighboring pixels, divided by the height of the center pixel. A pixel was classified as an “edge” if h > 0.5 m but was above the height CV threshold. The next steps were to merge all neighboring seed pixels into crown objects, grow each crown object by one pixel but only into edge pixels (to prevent merging of trees), and split oblong crowns with roundness > 0.6. These threshold values were iteratively developed until the automatically generated object outlines matched field-verified crown outlines visually delineated in the field. The final crown objects were labeled as “tree crowns” if their maximum H_obj was > 1.5 m or otherwise as “shrubs”. The projected crown area (m²) and maximum height (m) were calculated for each tree crown in eCognition.

2.6. BRDF Correction Model

The BRDF correction model used here was adapted to suit airborne hyperspectral data from the BRDF model initially developed for the MODIS reflectance/albedo product, the “Ross Thick, Li Sparse reciprocal combination”, hereafter referred to as the Ross-Li model [16,24,36]. Like the MODIS Ross-Li model, we used a semi-empirical, kernel-based model derived from radiative transfer approximations within a vegetation canopy. Our model is similarly composed of two kernels, a geometric scattering kernel (F₁) and a volume scattering kernel (F₂); however, the drastically different spatial resolution of airborne data compared to MODIS lead to the primary deviation from the Ross-Li model: the use of the Li-dense reciprocal kernel to compute F₁ (rather than Li-sparse). Preliminary testing of all four combinations of Ross-thin, Ross-thick, Li-sparse, and Li-dense showed that Li-dense with Ross-thick best reproduced target nadir spectra when given off-nadir spectra as input. This is supported by the physical consideration of individual leaves and branches within a 1.12 m pixel, which are likely to mutually shadow each other, thereby warranting the assumptions made in the Li-dense kernel. Despite the low leaf area index (LAI) in savannas relative to forests, Ross-thick, which assumes LAI >> 1, always predicted more accurate spectra than Ross-thin (LAI << 1). Thus, the modeled surface reflectance R was computed for each band as:

$$R\left({\theta }_s,{\theta }_v,\phi \right)=c_0+c_1{F}_1\left({\theta }_s,\hspace{0.17em}{\theta }_v,\hspace{0.17em}\phi \right)+c_2{F}_2\left({\theta }_s,\hspace{0.17em}{\theta }_v,\hspace{0.17em}\phi \right)$$

(1)

where θ_s is the solar zenith angle (radians), θ_v is the view zenith angle (radians), φ is the relative azimuth angle (radians), F₁ is the Li-dense kernel, F₂ is the Ross-thick kernel, and c_i are constants inverted from measurements, discussed below (where c₀ represents the isotropic reflectance). The kernels were computed from [24,37] (with the reciprocal modification to the Li kernel) as:

$$\begin{array}{c}{F}_1=\frac{\left(1+\text{cos}\hspace{0.17em}{\xi }^{\prime }\right)\hspace{0.17em}\text{sec}\hspace{0.17em}{\theta }_v^{\prime }\hspace{0.17em}\text{sec}\hspace{0.17em}{\theta }_s^{\prime }}{\text{sec}\hspace{0.17em}{\theta }_v^{\prime }+\text{sec}\hspace{0.17em}{\theta }_s^{\prime }-V}-2\\ V=\frac{1}{\pi }\left(t-\text{sin}\hspace{0.17em}t\hspace{0.17em}\text{cos}\hspace{0.17em}t\right)\hspace{0.17em}\left(\text{sec}\hspace{0.17em}{\theta }_v^{\prime }+\text{sec}\hspace{0.17em}{\theta }_s^{\prime }\right)\\ \text{cos}\hspace{0.17em}t=\frac{h}{b}\hspace{0.17em}\frac{\sqrt{D^2+(\text{tan}\hspace{0.17em}{\theta }_s^{\prime }\hspace{0.17em}\text{tan}\hspace{0.17em}{\theta }_v^{\prime }\hspace{0.17em}\text{sin}\hspace{0.17em}\phi }}{\text{sec}\hspace{0.17em}{\theta }_v^{\prime }+\text{sec}\hspace{0.17em}{\theta }_s^{\prime }}\\ D=\sqrt{{\text{tan}}^2\hspace{0.17em}{\theta }_s^{\prime }+{\text{tan}}^2\hspace{0.17em}{\theta }_v^{\prime }-2\hspace{0.17em}\text{tan}\hspace{0.17em}{\theta }_s^{\prime }\hspace{0.17em}\text{tan}\hspace{0.17em}{\theta }_v^{\prime }\hspace{0.17em}\text{cos}\hspace{0.17em}\phi }\\ \text{cos}\hspace{0.17em}{\xi }^{\prime }=\text{cos}\hspace{0.17em}{\theta }_s^{\prime }\hspace{0.17em}\text{cos}\hspace{0.17em}{\theta }_v^{\prime }+\text{sin}\hspace{0.17em}{\theta }_s^{\prime }\hspace{0.17em}\text{sin}\hspace{0.17em}{\theta }_v^{\prime }\hspace{0.17em}\text{cos}\hspace{0.17em}\phi \\ {\theta }^{\prime }={\text{tan}}^{-1}\left(\frac{b}{r}\text{tan}\hspace{0.17em}\theta \right)\\ F_2=\frac{4}{3\pi }\hspace{0.17em}\frac{1}{\text{cos}\hspace{0.17em}{\theta }_s+\text{cos}\hspace{0.17em}{\theta }_v}\hspace{0.17em}\left[\left(\frac{\pi }{2}-\xi \right)\hspace{0.17em}\text{cos}\hspace{0.17em}\xi +\text{sin}\hspace{0.17em}\xi \right]-\frac{1}{3}\end{array}$$

(2)

The values for the parameters h/b and b/r (object shape and height) set to 2.0 and 10.0 after iterative testing demonstrated these values minimized the residual difference between test and predicted spectra. A multiplicative model was used to compute the final BRDF-corrected reflectance, ρ_c, by multiplying the observed reflectance, ρ, by the ratio of the modeled reflectance at the reference geometry (i.e., nadir view and roughly the mean solar zenith angle across all flight lines of 40°) to the modeled reflectance at the observed geometry:

$${\rho }_c=\rho \frac{R\left({\theta }_s=40°,\hspace{0.17em}{\theta }_v=0,\hspace{0.17em}\phi =0\right)}{R\left({\theta }_{s,\mathit{obs}},\hspace{0.17em}{\theta }_{v,\mathit{obs}},{\phi }_{\mathit{obs}}\right)}$$

(3)

The BRDF model (Equation (1)) constants c₁–c₃ were calibrated for each reflectance mosaic using several thousand samples randomly selected throughout the mosaic. The reflectance mosaic was first classified into the following light scattering classes (based on NDVI and NIR filters): well-lit vegetation, shaded veg, grass, soil, other (e.g., water), and clouds. A separate set of the three BRDF model constants were then calibrated for each scatter class and waveband. To test the hypothesis that minimizing BRDF effects would dampen intra- and inter-flight line reflectance artifacts in the species predictions, only spectral data (no LiDAR) were used to train and compare SVM models with and without BRDF-corrected input.

2.7. SVM Model Testing, Calibration, and Validation

We chose support vector machines (SVM) as the classification method for our analysis due to their superior performance in comparison with other classification methods (e.g., maximum likelihood, neural networks) [11,38]. SVM is a non-parametric classifier which has gained popularity in the machine learning and remote sensing communities [39] due to its ability to efficiently handle large input space and its relative insensitivity to the Hughes phenomenon [12]. This method divides classes in hyperspace using hyperplanes fit to the training dataset. Because classes are rarely linearly separable in the original feature space, the data are transformed into a kernel feature space of higher dimensionality. We used the radial basis function (RBF) kernel due to its good performance compared to other kernel functions [11] and the low number of parameters which must be fit. The equation for the RBF kernel is as follows:

$$K\left(x_i,\hspace{0.17em}{x}_j\right)=\mathit{exp}\left(-\frac{{\Vert x_i-x_j\Vert }^2}{2{\sigma }^2}\right)$$

(4)

where K(x_i,x_j) is the distance between points x_i and x_j in feature space, and σ is the kernel width parameter which controls the tradeoff between over- and under-fitting. A second parameter, C, must also be optimized which controls the tradeoff between complexity of decision rule and frequency of training error [40]. These parameters are fit through an exhaustive grid search; the fitted values of the parameters are those which produce the best cross-validated classification accuracy. All SVM calculations were performed using the “e1071” package [41] of the R statistical software [42].

Five SVM models were tested, two single-level SVM models and three stacked (two-level) SVM models, in order to separately quantify effects from including height and aggregating species predictions at the crown level. Stacked generalization [43], or “stacking”, is a method of classification wherein the output of one classification is used as the input for a second classification. Here we present a modified stacked classification method, in which we aggregate the results of a hyperspectral-based, pixel-level species classification with LiDAR-derived, crown-level structural information (maximum height and projected area).

Each model was run 100 times and the overall accuracy of crown predictions was estimated each time by five-fold crown-level cross-validation. In the validation procedure, one fifth of the crowns were removed from the training data to make the test crowns, the model was trained using the remaining four fifths of the crowns, and predictions were made for the test crowns. This is repeated five times such that all crowns were predicted using 80% of the crowns. The overall accuracy is the percentage of crowns correctly classified.

Single-level SVM models:

• S: An SVM was performed on image pixels using spectral data only. Each pixel was assigned a probability of belonging to each of the classes. Crown identities were predicted by averaging class probabilities over all pixels in a crown and taking the class with the highest probability (n = 14,998 pixels and 72 spectral bands as input variables).
• S + H: An SVM was performed on image pixels using spectral and LiDAR height data. This is the same as Model 1, but with the additional variable of pixel height. Crown identities were predicted in the same way as in Model 1.

Stacked SVM models:

• S: An SVM was first performed on image pixels using only spectral data, as in Model 1. The average probabilities of all pixels in a crown were then used as the input variables to a second SVM (n = 729 crowns and 15 class probabilities as input variables to the second SVM).
• S + MH: The same procedure as in Model 3, with the additional variable of maximum crown height added to the crown-level SVM.
• S + MH + A (final model): The same procedure as in Model 4 with the additional variable of crown area in the crown-level SVM. Crown area was defined as the number of pixels belonging to a tree crown, before filtering based on NDVI and mean NIR reflectance.

In constructing this set of models, we were interested in a few key issues. First, we sought a way to integrate information contained in the multiple pixels belonging to a single crown. Preliminary analysis showed that using the mean spectra over all pixels in a crown as input into a crown-level SVM had poor performance (approximately 54% overall accuracy)—much worse than predicting the pixels individually and then averaging the results for each crown. We also found that using the averaged class probabilities of a pixel-level SVM as input into a crown-level SVM performed well as a way of combining data at the crown-level. We therefore included Model 3 in the set for comparison with Model 1. Second, we wanted to determine the best way to incorporate the LiDAR height data. We included height data in two different models: as a pixel-level and a crown-level variable (Model 2 and Model 4, respectively). At the crown level, we chose maximum height over average height because it less dependent on canopy shape, and preliminary analysis showed less variance in maximum height than mean height among our species. Lastly, we compared these models with the “final model” (Model 5), which incorporated pixel-level information from the spectral data as well as the crown-level variables of maximum height and crown area.

3. Results

3.1. SVM Classification of Woody Plant Species

We found that using spectral data in a stacked SVM (Model 3) yielded higher accuracy than when crowns were predicted using the pixel-level results from a single SVM (Model 1); see Figure 3. This indicates that no accuracy is lost when combining pixel-level information at the crown-level in this way; in fact, accuracy increased significantly (P < 0.001, paired t-test). The LiDAR height data were found to increase model accuracy whether used as a pixel-level or a crown-level variable. However, the gain in accuracy compared to an equivalent model using spectral data only was greater when adding maximum crown height as a crown-level variable (a 3.3% increase in accuracy from Model 3 to Model 4) than when adding pixel height to the pixel-level SVM (a 1.4% increase in accuracy from Model 1 to Model 2). The full model (Model 5) had the highest average overall accuracy (76.7%, 15 species classes) of all the models, which incorporated spectral data, maximum crown height, and crown area. It performed significantly better than the next best model (Model 4) (P < 0.001, paired t-test), which did not include crown area. Comparison of the producer and user accuracies for each species is shown in Figure 4.

3.2. Spatial Performance of Automated Crown Segmentation and Species Predictions

The tree crown polygons derived from the automated crown segmentation of the LiDAR data were found to be well registered to the corresponding crowns in the hyperspectral data (Figure 5). The LiDAR data better resolved small trees and delineated overlapping crowns than the hyperspectral data. Although the smallest detectable crown in the LiDAR data was 2.2 m in diameter (two pixels wide), small differences in registration between the hyperspectral and LiDAR data often resulted in an insufficient number of samples located in crowns of this size. The threshold of at least three hyperspectral samples mitigated this issue but also resulted in some crowns larger than 2 m in diameter not being predicted. This can be seen by comparing Figure 5(b,c), where some of the smaller crowns with low sample sizes are predicted as “other”. Conversely, many samples passed the NDVI and NIR filters for vegetation and appeared pink in the false-color infrared in Figure 5(b), yet they were (correctly) not classified since their maximum height was below the 1.5 m threshold and thus they were not segmented as a tree crown from the LiDAR data.

3.3. Spectral Performance of BRDF Correction

Variation in vegetation reflectance was observed to have a substantial BRDF component, with NIR reflectance of a given target tree varying by as much as 40% (relative change) between minimum and maximum VZA of −19° to 23° (Figure 6). VIS reflectance exhibited slightly higher sensitivity to BRDF effects, increasing by up to 50% over the same VZA range. After applying the BRDF model to adjust reflectance to nadir view angle (VZA = 0°) and 40° solar zenith angle, the difference in reflectance between minimum and maximum VZA dropped to less than 5% in both the VIS and NIR regions. The BRDF model also attenuated a cross-track artifact at the 940 nm water absorption feature, introduced during the atmospheric correction step.

3.4. Effects of BRDF Correction on Species Prediction

The cross-track brightness gradients discernible in the original hyperspectral color infrared image (Figure 7(a)) were almost entirely absent in the BRDF-corrected data (Figure 7(b)) for any bands above 450 nm (low signal-to-noise ratios prevented the eight bands below 450 nm from being used). Effects of BRDF on species prediction were evident when the predicted probability map of each species using the uncorrected hyperspectral data (Figure 7(c,e,g)) was compared to maps using the same model trained on BRDF-corrected spectra and then applied to the corrected spectra (Figure 7(d,f,h)). Of the twelve species classes in this test, ten had flight line artifacts detectable in the probability maps derived from the uncorrected spectral data (examples shown in Figure 7(c,e,g)).

Prior to BRDF correction the probability of a given species covaried with the across-track brightness gradients of each flight line of the original image. Two of the species with artifacts were not abundant in this spatial subset and had a very low mean probability (i.e., low signal to noise ratio), yet the remaining eight species, including common species with high mean probabilities, exhibited twofold or more variation in probabilities across-track. With BRDF correction all twelve species classes had reduced variation, with eight species absent of flight line artifacts entirely (including three examples shown in Figure 7(d,f,g)). Four species classes still had artifacts but were less severe than prior to BRDF correction, including A.tor, C.imb, E.div, and “other”.

4. Discussion

Mapping species at large scales using imaging spectroscopy requires a classifier with good generalization that can also handle the large dimensionality of this spectral data type [40,44]. The SVM is widely considered to meet these requirements in that it is reasonably resistant to over-fitting when given a large number of input variables [12]. The results reported here further demonstrate the robustness of the SVM approach with a high final overall accuracy (76%) despite the large dimensionality of the dataset and the low LAI of savanna trees. We note that preliminary testing indicated that the single largest increase in accuracy (from 65% to 71% for hyperspectral only) was due to increasing the training/validation sample size from 290 to 729 crowns. This suggests the SVM classifier remained sensitive to the Hughes effect, and we cannot rule out that further increasing the target sample size above 25 samples/species might increase overall accuracy.

More broadly, we found the SVM classifier well suited as a second-level aggregator of crown maximum height and species probabilities, which improved accuracy between major size classes (e.g., shrub, dwarf tree, large tree). A second SVM allowed inclusion of tree maximum height rather than pixel-level height, providing an increase in accuracy from 73.8% to 76.5%, notable given that no additional inputs were required. However, given the additional cost and complexity of having both LiDAR and hyperspectral data for a given study area, we note the relative improvement over hyperspectral data alone was relatively minor (from 73.0% to 76.5%). Hyperspectral data enable species detection even in small, low LAI savanna trees because reflectance spectra capture the underlying differences in leaf chemistry and structure between species. Species classification with LiDAR alone is not feasible since many species have similar height ranges. Yet preliminary analysis showed crown segmentation using the LiDAR data was qualitatively more reliable, than spectral-based segmentation, particularly for small trees, and the final species maps presented here relied on LiDAR-based crown segmentation. LiDAR data also play an important secondary role in airborne taxonomic mapping by accurately measuring the 3-D location of each hyperspectral sample, which aids georectification and mosaicking of multiple flight lines of hyperspectral data. The precise viewing and illumination angles used during the BRDF correction of the hyperspectral data were also derived from LiDAR data, although this could be done without LiDAR.

Airborne imaging spectroscopy studies vary in their approach to addressing BRDF effects [26–29,45,46], although they are often corrected using an empirical model (e.g., the modified Walthall model [46]). However, empirical models often perform poorly in the radiometrically heterogeneous nature of savannas (e.g., soil, grass, shrubs, trees). Another common practice is to not apply BRDF correction at all but to analyze individual flight lines rather than a mosaic, which prevents revealing the contrast along flight line boundaries but actually exacerbates BRDF effects by leaving flight line edges exposed. We found that prior to BRDF correction, the NIR reflectance varied by ∼40% as a function of view zenith angle, but after correction NIR varied by less than 5% between the minimum and maximum observed SVA (Figure 6).

Comparison of the total accuracies for all species before and after BRDF correction (using spectral data only to prevent confounding effects from LiDAR) indicated nearly equivalent overall predictive performance (71.0% and 70.3%, respectively). Yet this general metric is not reflective of the increased accuracy of many individual species due to dilution from three species classes: two of the most abundant species (C.api and A.nig) accounted for over 20% of the sample size, yet they did not have flight line artifacts prior to BRDF correction and thus had no accuracy improvement from BRDF correction. Given that the “other” species class is an amalgamation of many species, its accuracy had a negligible improvement from BRDF correction, as expected, yet accounted for 27% of the sample size with no change in overall accuracy. Due to these three classes, roughly half of the BRDF test sample size did not initially have flight line artifacts, which diluted the contribution of improved predictions of individual species to the overall prediction accuracy. For these reasons we conclude the post-BRDF visual patterns observed in the species probability maps suggests the BRDF correction is more broadly beneficial during large-scale prediction than that reflected by the modest increase in overall accuracy.

5. Conclusions

We focused on the challenges facing ecosystem-scale species mapping using a combination of LiDAR and imaging spectrometer data. By automatically segmenting crowns and then applying a stacked SVM model, we addressed the issue of incorporating variables best expressed at different scales (e.g., pixel vs. crown) and from different sensors. The varying anisotropic properties found in savanna landscapes is recognized and modeled using a semi-empirical BRDF model, adapted from the MODIS BRDF model for application to high resolution airborne hyperspectral data. Our results demonstrate the high prediction accuracy that can be achieved using a SVM classifier, and that application of a BRDF model to airborne spectral data can greatly reduce flight line artifacts, despite the spectral heterogeneity of savanna landscape. We find that confronting these three challenges—anisotropy, integration of multi-scale information, and crown segmentation—enables species mapping at ecosystem scales, a critical starting point for monitoring biodiversity and ecosystem function.