Classification of Grassland Successional Stages Using Airborne Hyperspectral Imagery

2.1. Study Area

The study area (centred on 56°40′49″N, 16°33′58″E) covers approximately 22.5 km² and is located on the Baltic island of Öland in SE Sweden (Figure 1A). The bedrock consists of Cambro-Silurian limestone, the average elevation is approximately 36 m above sea level, and the overall topography is flat [48]. The area is crossed by a few low ridges of glaciofluvial deposits. The mean annual temperature is 7 °C and the mean annual precipitation is 468 mm [48]. The present-day landscape consists of a mosaic of arable fields, villages, forests, and grasslands. The majority of the grasslands are grazed, with varying intensity, mainly by cattle.

Detailed information on the historical land use within the study area is available [49]. Whereas the study area was dominated by grasslands during the 18th century [49], most (>80%) of the ancient grasslands have been transformed to arable land during the last 300 years [49]. Approximately 15% of the current grasslands have developed during the last 50 years on previously arable land [49].

2.2. Grassland Sites and Grassland Age-Classes

Using historical and present-day land-use maps, aerial photographs [49] and field inventories, we identified 299 grazed grassland sites that were separated from each other and the surrounding landscape by walls or fences. With the help of geographic information system (GIS) overlay analysis of the land-use maps and aerial photographs, the sites were categorized according to their grazing continuity (grassland age) and assigned to one of three age-classes within the arable-to-grassland succession; 5–15 (young grasslands), 16–50 (intermediate-aged grasslands), and >50 (old grasslands) years of grassland continuity (Figure A1).

Within each site, two sampling points were randomly positioned in open (not covered by shrubs or trees) grassland vegetation, with the constraints that they should be at least 25 m apart (to minimize effects of spatial autocorrelation in the vegetation), at least 13.5 m from the site boundary (to minimize edge-effects in the vegetation [50]), and at least 13.5 m from shrubs or trees higher than 50 cm (to minimize shading-effects). In total 239 grassland sites (89 young, 73 intermediate, and 77 old) out of the 299 sites could accommodate these constraints. From these 239 grassland sites, we randomly selected 60 sites (20 young, 20 intermediate, and 20 old). Within these 60 sites, a bioassay approach (cf. [51,52]) based on indicator species (such as Sesleria caerulea and Molinia caerulea) was used to define sites with “dry” grassland vegetation, and exclude moister grassland vegetation. A total of 52 sites (17 young, 18 intermediate, and 17 old) out of the 60 sites were characterized by dry grassland vegetation and used for the vegetation and remote sensing sampling (Figure 1B). Of the 17 “old” sites, 13 sites had a management continuity of >280 years [49]: the vegetation in these sites falls within the Natura 2000 habitat type “Fennoscandian lowland species-rich dry to mesic grasslands”.

A hand-held differential global positioning system (DGPS) receiver (Topcon GRS-1 GNSS, equipped with a PG-A1 external antenna (Topcon Corporation, Japan)) connected to a real-time positioning service (SWEPOS) was used to log (with an accuracy of ∼1 cm) the ground coordinates of the two sampling points within each of the 52 grassland sites. The sampling points were divided into two data sets (a training and a validation data set) by randomly assigning the two sampling points from each site to one or other of the two data sets (Figure 1C).

2.3. Plant Community Characteristics

Vegetation sampling was carried out between 15 May and 15 July 2011. A 1 m × 1 m plot was centred over each of the two sampling points within each of the 52 sites. Each of the 104 (1 m × 1 m) plots was divided into a grid of 100 (10 cm × 10 cm) sub-plots, and the presence of all non-woody vascular plants was recorded within each sub-plot (Table A1). For each 1 m × 1 m plot, the frequency of each species was calculated as the number of sub-plots in which the species was present.

The plant species recorded in the 1 m × 1 m plots were assigned values for SLA and two Ellenberg indicator values (for nutrients and moisture availability). The trait information was compiled from the LEDA trait database [53], which provides species-level SLA values that represent aggregated trait values. The LEDA values are based on multiple individuals per species, measured under different environmental conditions. Trait data were not available for the full set of species (n = 185) in the present study. For SLA, previous grassland studies have shown that imputation methods based on multiple imputation by chained equations may be used for filling gaps in functional trait databases [54]. Following Taugourdeau et al. [54], estimates for missing values (approximately 9% of the species) were obtained using the multivariate imputation by chained equation (MICE) method [33,55,56]. Ellenberg indicator values were extracted from the JUICE database [57]: values for the four species that were not found in the database were extracted from Ellenberg et al. [58]. A frequency-weighted mean value (CWM) [35] was calculated, for SLA and each of the Ellenberg values, for each plant community: CWM(x) = ∑_ip_i × x_i, where p_i is the relative frequency of the ith species and x_i is the trait or Ellenberg index value of the iith species.

The percentage ground cover of grasses was estimated within four 50 cm × 50 cm sub-plots within each 1 m × 1 m plot, and a mean within-plot value for grass cover (excluding dead litter) was calculated for each plot.

2.4. Remote Sensing Data

Remote sensing data over the study area were acquired on 9 July 2011. Twenty-five flight lines were recorded at around solar noon, using two airborne HySpex hyperspectral spectrometers (Norsk Elektro Optikk AS (NEO), Lörenskog, Norway), and a push-broom scanning mode at a flight altitude of approximately 1500 m. The flight was carried out by Terratec AS, Lysaker, Norway. All flight lines were conducted either from north to south or from south to north to minimize illumination effects. The output from the two HySpex spectrometers (VNIR-1600 operating over the 414 to 991 nm range and SWIR-320m-e operating over the 966 to 2501 nm range) consisted of 416 wavebands with spectral resolutions of 3.7 nm (VNIR-1600) and 6.0 nm (SWIR-320m-e) (Figure 2). The spatial resolution (pixel size) of the image captured from the VNIR-1600 spectrometer was 0.5 m × 0.5 m and the spatial resolution of the SWIR-320m-e image was 1.0 m × 1.0 m. A digital elevation model (DEM) was created with the help of LiDAR data recorded during the flight.

2.5. Statistical Analyses

3.1. Classification Based on the Full Set of 269 HySpex Wavebands (Model 1)

The misclassification rate decreased progressively with the number of latent variables (from 57.7% with one latent variable), reaching a minimum (30.8%; Figure 3) with the first four latent variables before increasing again. We used the first four latent variables in the final model.

The producer’s and user’s classification accuracies for the individual grassland age-classes varied between 47% and 100% for the training data set, and between 65% and 94% for the validation data set (Table 1).

The overall classification accuracy and the Kappa statistic value were 77% and 0.65, respectively, for both the training and the validation data set (Table 1).

3.2. The HySpex Wavebands That are Most Important for the Classification of Grassland Age-Classes

The VIP value is an indicator of the relative influence of each predictor in a PLS model. Out of the 269 HySpex bands used as explanatory variables in the PLS-DA model, 177 bands had VIP values greater than 0.8, indicating that these bands were the most important predictor variables in the remote sensing-based classification of grassland age-classes (Figure 4). In total, 14 wavebands in the blue region (414–499 nm), 16 wavebands in the green region (505–596 nm), 20 wavebands in the red region (602–716 nm), 17 wavebands in the red-edge portion (722–818 nm), 23 wavebands in the NIR1 part (824–969 nm), 39 wavebands in NIR2 part (975–1394 nm), and 48 wavebands in the SWIR portion (1448–2417 nm) of the electromagnetic spectrum were associated with VIP scores greater than 0.8 (Figure 4).

3.3. Classification Based on a Reduced Set of HySpex Wavebands (Model 2)

The subset of HySpex bands used to develop Model 2 included the 177 spectral bands that were associated with VIP values greater than 0.8 in the model based on the full set of 269 HySpex bands. The inclusion of four latent variables gave the lowest misclassification rate in the second PLS-DA model.

The producer’s and user’s classification accuracies for the individual grassland age-classes for Model 2 varied between 65% and 94% for the training data set (Table 2). The classification accuracies for the individual grassland age-classes for the validation data set ranged between 74% and 94% (Table 2). The overall classification accuracy and the Kappa statistic value were 81% and 0.71, respectively, for the training data set and 85% and 0.77, respectively, for the validation data set (Table 2).

For the validation data sets, the overall classification accuracy and the Kappa statistic value were higher for Model 2 than for Model 1 (Tables 1 and 2). The differences in classification accuracies between the two models were 8% for the overall accuracy and 0.12 for the Kappa statistic value (Tables 1 and 2).

3.4. Plant Community and Spectral Characteristics of Grasslands Representing Different Age-Classes

The first four latent variables (LV1–LV4) in Model 2 explained 34%, 41%, 12%, and 10% of the variation in the explanatory matrix (spectral data), respectively. Old (negative scores) and intermediate-aged (positive scores) grassland plots were separated on the first latent variable (LV1), whereas LV2 discriminated between young (positive scores) and intermediate-aged (negative scores) plots (Figure 5A).

Correlation coefficients between the scores for grassland plots on the first two latent variables, and the ecological variables (SLA, Ellenberg indicator value for nutrients and moisture availability, and the percentage cover of grasses) were calculated with the help of multiple regressions. Significance values for the correlation coefficients were based on 9999 permutations. The correlation coefficients are represented as arrows in Figure 5A, with the direction and length of the arrow indicating, respectively, the sign and the strength of the correlation coefficient. The scores for the grassland plots were significantly associated with the mean within-plot Ellenberg indicator values for nutrient availability (p = 0.010), and with the within-plot cover of grasses (p = 0.045). The young plots were characterized by the highest and the old plots by the lowest indicator values for nutrients, whereas the intermediate-aged plots were characterized by the highest cover of grasses (Figure 5A). The mean within-plot SLA and Ellenberg indicator value for moisture were not significantly associated with the plot scores in the PLS-DA.

The old plots were associated with negative loadings for the blue, red, and SWIR regions on the first latent variable, and the intermediate-aged plots were associated with positive loadings for the red-edge and NIR1 regions (Figure 5B). The old plots, which are characterized by lower indicator values for nutrients than the intermediate-aged plots (Figure 5A), had higher reflectance values than the intermediate-aged plots in the blue, red, and SWIR regions. The intermediate-aged plots, which are characterized by higher cover of grasses than the old plots (Figure 5A), had higher reflectance values than the old plots in the green, red-edge, and NIR1 regions (Figure 5B).

On LV2, the young plots were associated with positive loadings for the red-edge, NIR1, and NIR2 regions, while the intermediate-aged plots were associated with negative loading for the green region (Figure 5B). The young plots, which are characterized by higher indicator values for nutrients than the intermediate-aged plots (Figure 5A), had higher reflectance values than the intermediate-aged plots in the red-edge, NIR1, and NIR2 regions (Figure 5B). The intermediate-aged plots, which are characterized by higher cover of grasses than the young plots (Figure 5A), had higher reflectance than the young plots in the green region (Figure 5B).

Species lists for each of the three grassland age-classes are provided in Table A1 (Supporting information). Figure A1 (Supporting information) shows a field photograph of a sample plot representing each of the three grassland age-classes.

4.1. PLS-DA Based Classification of Grassland Age-Classes

The three age-classes in our study were successfully classified using a pre-selected set of hyperspectral wavebands and PLS-DA, indicating that fine-scale hyperspectral measurements are able to accurately identify grassland successional stages within a local landscape. Our results are consistent with those of earlier studies of other types of vegetation that also show that the exclusion of wavebands that provide little information related to the response variable improves the PLS-DA based classification of plant communities [47,81]. For example, using airborne AISA Eagle hyperspectral imagery, Peerbhay et al. [47] showed that a PLS-DA model based on the 78 wavebands that were most relevant for the classification of forest species gave 8.17% higher overall classification accuracy than a model utilising all 230 AISA Eagle wavebands. Our results are also consistent with the results of Peerbhay et al. [47] in that they show that a high number of hyperspectral wavebands can be compressed into a few latent variables with the help of PLS, reducing the risk of model overfitting and, at the same time providing a successful classification of different vegetation types.

4.2. Spectral Dissimilarities between Grasslands Belonging to Different Age-Classes

The majority of the present-day grazed grasslands in Sweden originate from abandoned arable fields that have been cultivated as leys [82] and seeded with grasses, such as Lolium perenne, Festuca pratensis, Poa pratensis, and Phleum pratense (cf. [83]). The effects of the addition of nutrients during arable cultivation may persist during the early stages of the arable-to-grassland succession [84,85]. Recently cultivated fields, such as the young grasslands in our study (Figure 5A), are generally characterized by a high availability of nutrients that allows nitrophilous and ruderal species (fast growing and disturbance-tolerant species) to invade and restricts the establishment of less competitive species during early succession (cf. [33]). Early-successional grasslands are typically characterized by species that are adapted to nutrient-rich environments (e.g., Dactylis glomerata, Poa pratensis, and Geranium molle [86]) and nutrient-rich vegetation commonly contains plants with a relatively high SLA [87]. An increased availability of nutrients may result in greater above-ground production of biomass (cf. [85]), which, in its turn, influences the light conditions within the vegetation canopy: an increased leaf area represents a common adaptation to low light environments [38,88].

Residual effects of fertilization during arable cultivation may still persist within the intermediate-aged grasslands in our study (cf. [84,85]). However, the fact that continuous grazing management contributes to the removal of nutrients explains the reduced abundance of nitrophilous species during the second successional time step. Previous investigations within the study area have shown that grasslands with a grazing continuity of 16–50 years (corresponding to the intermediate-aged grasslands in the present study) are characterized by significantly lower community mean values for SLA than grasslands that have been grazed for 5–15 years (corresponding to the young grasslands in the present study) [33].

During the third successional time step, a long continuity of grazing management has contributed to a substantial reduction of nutrient levels, allowing species with a low competitive ability to establish in the old grasslands. Purschke et al. [33] showed that, within the study area, values of community-level SLA are lower in grassland sites that have been grazed for more than 51 years (corresponding to the old grasslands in the present study) than in sites with a shorter grazing continuity. Long grazing continuity is essential for the maintenance and establishment of species-rich ancient grassland vegetation characterized by, for example, Festuca ovina and Helianthemum nummularium [89].

Changes in habitat conditions during the arable-to-grassland succession have been shown to be accompanied by changes in plant functional composition and vegetation properties within the same study area [33,90]. In addition, changes in, for example, leaf nitrogen content, SLA, leaf chlorophyll content, leaf water content, and above-ground biomass have shown to be followed by changes in the vegetation spectral response within a range of different habitats [26,91–93]. For example, previous remote sensing studies have revealed that differences in leaf nitrogen content have strong influence on the reflectance from vegetation across the electromagnetic spectrum (visible (400–700 nm) [94], NIR (700–1300 nm) [95], SWIR (1300–2500 nm) [96]). The red-edge region has, been used to estimate the nitrogen concentration in ryegrass (Lolium spp.) [97], and SLA has previously been shown to influence reflectance in several wavebands. For example, Lymburner et al. [98] estimated SLA at the landscape scale with the help of broadband satellite data, and Asner et al. [99] showed that hyperspectral measurements were related to the SLA of tropical forest leaves. The SWIR wavelengths have been shown to be most important for estimation of leaf mass per area (the inverse of SLA [100]).

In the present study, differences in ecological variables (availability of nutrients and cover of grasses) between grassland age-classes are associated with differences in the spectral response from the grasslands (Figure 5). The low blue (414–499 nm) as well as the red (602–716 nm) and the high red-edge (722–818 nm) reflectance for the young plots (Figure 5B) suggest that these plots have a higher chlorophyll content in the vegetation than the old plots. Earlier remote sensing studies have shown that the chlorophyll content in grasslands can be assessed at both fine and landscape scales [46,101], and Darvishzadeh et al. [46] used the red-edge region to map canopy chlorophyll content in heterogeneous grasslands.

The higher cover of grasses (excluding dead litter) in the intermediate-aged plots compared with the young and old plots can be expected to lead to higher green (505–596 nm) reflectance in the intermediate-aged plots. Spectral measurements near the green reflectance peak at 550 nm have previously been shown to be important for the classification of vegetation communities (e.g., [47]).

The high red-edge (722–818 nm) and NIR1 (824–969 nm) reflectance for the young and intermediate-aged plots suggest that these plots were associated with larger above-ground biomass than the old plots (cf. [45]). The low NIR2 (975–1394 nm) reflectance for the intermediate-aged and the low SWIR (1448–2417 nm) reflectance for the young plots indicate that these were characterized by higher canopy moisture than the old plots. The SWIR wavelengths are strongly affected by water absorption [100]; in particular, the moisture-sensitive bands around 1200 nm (located in the NIR2 region in the present study) have been shown to be sensitive to the water content of vegetation in mixed woodland [102] and in grasslands [68].

Trampling by grazing animals is likely to lead to differences in the cover of bare soil between the different age-classes. Spectral dissimilarities between the grasslands belonging to different successional stages may also be explained by between-site edaphic dissimilarity (cf. [103]). The red part of the visible electromagnetic spectrum has been shown to be highly influenced by variation in soil properties [104]. Furthermore, increased levels of shading are expected within the vegetation cover as the above-ground biomass increases. Increased amounts of shadow are accompanied by a decrease in reflectance in all wavebands [105]. The presence of shadows may also have contributed to between-stage spectral dissimilarities in the present study. Monocotyledonous plants, such as grasses, have a compact leaf mesophyll structure, and a lower reflectance in the NIR than plants with a more porous mesophyll structure [106]. The high cover of grasses within the intermediate-aged grasslands in our study may have contributed to between-grassland stage spectral differences in the NIR. Grazing management decreases the accumulation of dead above-ground biomass [107] and the amount of litter has shown to affect the spectral response of grassland vegetation [108]. Our different-aged grassland plots—characterized by different continuity of grazing—may be expected to be associated with different amount of litter, which, in its turn, may have contributed to the spectral dissimilarities between the plots.

The results of our study are consistent with the results from previous ecological studies, which show that plant community characteristics such as the cover of grasses and community-weighted mean Ellenberg indicator values for nutrient availability, change during grassland succession [34,109]. The results of the present study also agree with those of previous interdisciplinary studies within remote sensing and ecology. Earlier studies have shown, for example, that hyperspectral reflectance measurements in grasslands are related to community-weighted mean Ellenberg indicator values for nutrient availability and for plant functional traits, such as SLA [12,32,93].

Although the accuracy of the predictions we obtained using models 1 and 2 were high, some grassland plots in our sample study were not classified correctly by either model. There are a number of possible reasons for this discrepancy. First, we conducted the spectral measurements in early July, when leaf senescence in response to summer drought is likely to have reduced the spectral differences between grassland age-classes. Second, hyperspectral images are known to be affected by radiometric noise [110]. While we used the commonly applied Savitzky-Golay filter [63], other denoising methods, such as the Minimum Noise Fraction (MNF) transform [111], might have improved the noise-reduction in our data. Thirdly, PLS-DA may not be the most optimal method for classifying grassland age-classes. Other classification methods, such as artificial neural networks or support-vector machines (cf. [112]), which are able to take into account non-linear relationships between dependent and explanatory variables, might provide better discrimination between grassland successional stages.

4.3. Future Directions

Future research targeting the improvement of hyperspectral methods for classifying grassland successional stages may include the development of vegetation spectral indices (e.g., [93]), which could potentially provide enhanced information on the plant community characteristics associated with different grassland successional stages. There is also a need for studies that focus on the development of classification methods for discriminating between age-classes at the scale of entire grassland sites. In addition, in order to be able to develop the possible use of hyperspectral data as a source of information within ecological research, we need a better understanding of the ways in which different plant community variables (characterizing different grassland successional stages) influence the hyperspectral response. The spectral response of vegetation canopies changes over the vegetation period as a result of temporal variation in the biochemical and structural properties that influence the reflectance of the vegetation canopy [44]. Additional studies will need to explore seasonal variation in the relationships between hyperspectral data and plant community variables.