Detecting shrub encroachment in seminatural grasslands using UAS LiDAR.

Shrub encroachment in seminatural grasslands threatens local biodiversity unless management is applied to reduce shrub density. Dense vegetation of Cytisus scoparius homogenizes the landscape negatively affecting local plant diversity. Detecting structural change (e.g., biomass) is essential for assessing negative impacts of encroachment. Hence, exploring new monitoring tools to achieve this task is important for effectively capturing change and evaluating management activities.This study combines traditional field-based measurements with novel Light Detection and Ranging (LiDAR) observations from an Unmanned Aircraft System (UAS). We investigate the accuracy of mapping C. scoparius in three dimensions (3D) and of structural change metrics (i.e., biomass) derived from ultrahigh-density point cloud data (>1,000 pts/m2). Presence-absence of 12 shrub or tree genera was recorded across a 6.7 ha seminatural grassland area in Denmark. Furthermore, 10 individuals of C. scoparius were harvested for biomass measurements. With a UAS LiDAR system, we collected ultrahigh-density spatial data across the area in October 2017 (leaf-on) and April 2018 (leaf-off). We utilized a 3D point-based classification to distinguish shrub genera based on their structural appearance (i.e., density, light penetration, and surface roughness).From the identified C. scoparius individuals, we related different volume metrics (mean, max, and range) to measured biomass and quantified spatial variation in biomass change from 2017 to 2018. We obtained overall classification accuracies above 86% from point clouds of both seasons. Maximum volume explained 77.4% of the variation in biomass.The spatial patterns revealed landscape-scale variation in biomass change between autumn 2017 and spring 2018, with a notable decrease in some areas. Further studies are needed to disentangle the causes of the observed decrease, for example, recent winter grazing and/or frost events. Synthesis and applications: We present a workflow for processing ultrahigh-density spatial data obtained from a UAS LiDAR system to detect change in C. scoparius. We demonstrate that UAS LiDAR is a promising tool to map and monitor grassland shrub dynamics at the landscape scale with the accuracy needed for effective nature management. It is a new tool for standardized and nonbiased evaluation of management activities initiated to prevent shrub encroachment.

INTRODUCTION

Plant diversity in natural and seminatural grasslands experiences severe pressure from shrub encroachment (Timmermann, Damgaard, Strandberg, & Svenning, 2015; Wilsey, 2018). Establishment of woody species is a successional stage gradually occurring in grasslands and may induce encroachment if not exposed to natural processes, such as browsing, trampling and fires, or if grasslands are not managed otherwise. (D'Odorico, Okin, & Bestelmeyer, 2012). Shrub encroachment is enhanced by nutrient enrichment (Stevens, Dise, Mountford, & Gowing, 2004) and land‐use abandonment (Deák, Valkó, Török, & Tóthmérész, 2016) that leads to a decrease in plant species richness unless management strategies are implemented to keep shrub dominance low (Kesting, Petersen, & Isselstein, 2015). Grassland plant diversity is generally highest at intermediate levels of biomass (Ejrnæs & Bruun, 2000), and this relationship is also found for shrub encroachment (Kesting et al., 2015). High density of Cytisus scoparius is known to homogenize the landscape and to decrease plant diversity, especially in grasslands (Bellingham, 1998). The presence of C. scoparius can favor the establishment of fast‐growing species by decreasing the soil C:N ratio and hereby increasing N availability (Haubensak & Parker, 2004). Cytisus scoparius vegetation is therefore considered problematic to grasslands including many areas within its European native range (Sheppard, Hodge, Paynter, & Rees, 2002). The number of C. scoparius individuals in an area is not necessarily representative for invasiveness, and Parker (2000) suggests that structural changes in C. scoparius (e.g., biomass and density) are important to assess the negative impact on plant diversity. Hence, the structural components are necessary for evaluating the risk of grassland plant diversity decline resulting from encroachment and to assess the efficiency of management activities. Monitoring efforts that highlight the need for action (e.g., shrub reduction) should be strategically aligned with the actual spatial scale of management (Magurran, 2016). However, monitoring shrub dynamics with the detail (i.e., small grain size) and spatial extent needed for management is challenging (Cao, Liu, Cui, Chen, & Chen, 2018).

Remote sensing covers larger areas than classical field assessments of structural change, which are time‐consuming and provide only local information (Wachendorf, Fricke, & Möckel, 2017). Satellite‐based spectral information can provide valuable information on cover of C. scoparius at the landscape scale for areas with high shrub densities, although limited to observations in the flowering period (Hill, Prasad, & Leckie, 2016). However, the relatively low resolution of satellite imagery (Aplin, 2005) does not allow the quantification of structural information, for example, height and biomass, from individuals or lower concentrated areas of shrubs. In contrast, observations from cameras or LiDAR (Light Detection And Ranging) sensors mounted on an Unmanned Aircraft System (UAS, also referred to as UAV or drone) can provide very high‐resolution information on canopy (Wallace, Lucieer, Malenovský, Turner, & Vopěnka, 2016) or vegetation structure (Forsmoo, Anderson, Macleod, Wilkinson, & Brazier, 2018; Wang et al., 2017).

LiDAR technology is a remote sensing method providing three‐dimensional (3D) point cloud information suitable for quantifying vegetation structure (Lefsky, a, Cohen, W. B., Parker, G. G., & Harding, D. J., 2002). Airborne LiDAR has been used to map grassland vegetation (Zlinszky et al., 2014) and to explain variation in species diversity across varying vegetation communities (Moeslund et al., 2019). Especially, forestry research has successfully implemented airborne LiDAR data to, for example, detect the composition of gymnosperm species in plantations (Donoghue, Watt, Cox, & Wilson, 2007), to evaluate ecosystem services (Vauhkonen, 2018) and management efforts based on forest structure (Valbuena, Eerikäinen, Packalen, & Maltamo, 2016), or to demonstrate how LiDAR can be used to measure vegetation height of sagebrush (Mitchell et al., 2011). Furthermore, stages of shrub encroachment and biomass estimates have been mapped on a coarser resolution (30 m raster) based on a LiDAR point density of 5.6 points/m2 (Sankey, Shrestha, Sankey, Hardegree, & Strand, 2013). With advancing technology, it is now possible to mount good quality scanners on UAS (Manfreda et al., 2018). UAS LiDAR systems provide new opportunities to provide ultrahigh point density (>1,000 pts/m2) on demand; hence, mapping vegetation structure with high detail and sampling frequency becomes possible across areas of up to several square kilometers. UAS LiDAR systems have been used to detect individual trees and to measure metrics such as height and stem diameter (Wallace, Lucieer, Watson, & Turner, 2012; Wieser et al., 2017). The ability to detect single trees has been found to increase with higher point densities (Wallace, Lucieer, & Watson, 2014), and Balsi, Esposito, Fallavollita, & Nardinocchi, 2018 could contrast different shapes of horizontal overlapping trees using information from the whole volume of points in the point cloud. Furthermore, Moeslund et al., 2019 analyzed the potential of using airborne LiDAR‐derived metrics, including a biomass measure, to assess the diversity of different organisms (i.e., vascular plants, fungi, lichens, and bryophytes).

The aim of this study was to assess the value of UAS LiDAR for monitoring structural change in a seminatural grassland threatened by encroachment of C. scoparius. We developed a semiautomatic workflow to measure structural features of different shrubs and to enable standardized monitoring for estimating spatiotemporal changes of C scoparius biomass. Specifically, we addressed the following questions: (a) How accurately can C. scoparius be classified in UAS LiDAR point clouds? (b) How precise can the biomass of C. scoparius individuals be estimated from point clouds? (c) How does the estimated C. scoparius biomass change between two surveys at different seasons?

MATERIALS AND METHODS

Study area

The study was conducted in a seminatural grassland located within Nationalpark Mols Bjerge, Denmark (56°13′40.0″N; 10°34′30.2″E). The area is situated in a temperate climate zone with a mean annual temperature of 7.5°C and mean annual precipitation of 585 mm (Fick & Hijmans, 2017). The terrain elevation varies between 30 and 60 m above sea level. Graminoids and small, broad‐leaved herbs characteristic for European dry grasslands dominate the vegetation, while single standing trees and different species of shrubs are patchily distributed across the area. Cytisus scoparius forms dense stands in parts of the area and a broad‐leaf forest is located toward the northeast. Besides naturally occurring wildlife, the area is grazed by Galloway cattle and Exmoor ponies and thus affected by all year grazing and trampling. The grazers were introduced in 2016 as part of a rewilding experiment for creating a more self‐regulating ecosystem (Svenning et al., 2016).

UAS LiDAR system

We integrated the Surveyor laser‐scanning system (YellowScan, Montferrier sur Lez, France) on an 8‐rotary wing UAS (MK8‐3500; Mikrokopter, HiSystems GmbH; see Figure 1). The Surveyor is a dual‐return system ranging in the 903 nm wavelength with 360 degrees scanner angle, based on the Velodyne VLP‐16 “Puck” laser scanner with a maximum measurement range of 100 m and a ranging accuracy of 3 cm. The LiDAR sensor payload weighs 1.6 kg. The Surveyor has a global navigation satellite system (GNSS) receiver and inertial measurement unit (IMU) integrated that acts as a rover (Applanix APX15). We utilize a Trimble base–rover postprocessed kinematic (PPK) solution to gain subdecimeter accuracy for the system in XYZ directions (Chaponnière & Allouis, 2016).

FIGURE 1

UAS LiDAR system with background of Cytisus scoparius shrubs. Photo credit: Urs A. Treier

Data acquisition

Fieldwork was conducted in October 2017 (leaf‐on) and April 2018 (leaf‐off). Each time the UAS LiDAR system surveyed the area in three separate flight paths planned with the software Kopter‐Tool V2.20b with a flying height of 40 m above ground and a ground speed of 3 m/s. The planned waypoints ensured a standardized flight plan in transects with approximately 15 m between flight paths. We restricted the scan angle to ±55° (roll axis) with respect to the sensor pointing nadir (−90° pitch axis) resulting in an estimated overlap of 65% and thus reducing noise associated with very oblique scan angles.

Between the two flight campaigns, the position of 180 individuals of 12 different shrub/tree taxa (65 of these were C. scoparius) was determined with an RTK GNSS receiver on the ground (<2 cm absolute accuracy) for training and validation data. The additional genera measured, included Juniperus, Rubus, Rosa, Quercus, Betula, Pinus, Sambucus, Crataegus, Prunus, Malus, and Calluna, all collected along random transects with complementary samples of less abundant species outside (Appendix 2).

To validate the estimation of biomass from C. scoparius, ten individuals were harvested for dry weight measurements after 2 days of incubation at 60°C. Before cutting the shrubs, they were measured with the RTK GNSS device, resulting in small manual point clouds of 25 points each (Appendix 3).

Processing workflow

We developed a workflow for processing ultrahigh‐density point clouds from UAS LiDAR to detect and map structural change in shrubs (Figure 2). Point cloud data processing was performed separately for 2017 and 2018 data with OPALS software v. 2.3.1 specifically developed for handling airborne LiDAR data (Pfeifer, Mandlburger, Otepka, & Karel, 2014).

FIGURE 2

Flowchart of UAS LiDAR processing chain to derive structural information (biomass) from classified shrub species and comparing between two time periods

Preprocessing and quality control of LiDAR data

The position data recorded during the flight by the LiDAR system were postprocessed using PosPac UAV v. 8.2 with data recorded by a Trimble base station to obtain PPK corrected trajectory information and point clouds. All flight lines were carefully considered and selected using a QGIS (QGIS Development Team, 2019) plug‐in provided by YellowScan (YellowScan, 2016) to minimize noisy data from UAS turns. The resulting point clouds were then evaluated with 12 RTK GNSS measurements of fenceposts functioning as control points for the vertical accuracy (RMSE: 3.9 cm in 2017 and 6.1 cm in 2018) and, likewise, five ground markers in 2018 to evaluate horizontal accuracy (RMSE: 1.4 cm).

We performed a relative adjustment of the point cloud to improve the alignment of objects (i.e., vegetation). The procedure was based on fitting overlapping flight strips with a least squares matching technique as described in Ressl, Mandlburger, and Pfeifer (2009). This process aligned the point clouds better relatively by, for example, reducing the variation in ground point height. However, because of the variation in flight strips, the absolute accuracy decreased to a vertical RMSE of 7.6 cm and 9.1 cm for 2017 and 2018, respectively, and a horizontal RMSE in 2018 of 6.5 cm. For further details on preprocessing and quality control, see Appendix 1.

Applying thresholds to obtain a shrub layer point cloud

We coarsely divided the point cloud into height‐based vegetation classes and a ground class to reduce processing time for the later shrub classification. A digital terrain model (DTM) derived from the minimum point height within a 1x1 m moving window was used to calculate normalized height of points (NormZ) (Appendix 2). We then classified the points using thresholds into ground (NormZ < 0.15 m), low vegetation (NormZ < 0.3 m), and high vegetation (NormZ > 3.5 m) leaving a shrub layer point cloud of medium height vegetation (0.3 m < NormZ < 3.5 m). Calluna vulgaris fell into the low vegetation layer and therefore was excluded from the further shrub classification procedure.

Generating point‐based reference data of structural features

From the 180 measured GNSS points, we gathered adjacent LiDAR points within a size‐specific area to comprise a reference dataset of 13 classes (11 shrub taxa, fenceposts, and shrub absence points; Appendix 2). Reference data were collected as presence/absence (PA) of shrubs in late 2017 and early 2018, that is, between the UAS flights. Within this period, we would only expect shrubs to disappear, for example, browsing/trampling by animals, which would mean that such individuals will be included as a GNSS record in 2017 but would not be there in 2018. Because of the digital assignment of training/validation data, no points would represent the shrub in 2018, and therefore, no points will be assigned to the given shrub class.

We developed structural features derived from the UAS LiDAR point cloud to represent vegetation morphology with the aim to classify shrub species (Appendix 2). We targeted the variables to represent shrub species on varying levels of scale, ranging from small‐scale leaf characteristics to larger‐scale shrub growth form and shape. The shrub growth form is here interpreted as the general appearance of the shrub and of how branches and leaves are arranged. Hence, it corresponds to specific, although overlapping features, such as structural complexity, density, and light penetration (Popescu, 2011). Specifically, we calculated the variance in height from a fitted normal‐based plane (Pseudowaveform) to resemble roughness or structural complexity (Van Aardt et al., 2012). The amount of points (Point count) and average distance between points (Point distance) were used as a measure of density, while light penetration was represented by including the count of ground points classified by the above mentioned threshold (Ground points count) and the average number of returning echoes (No. of echoes). Additionally, we calculated a Rank feature, where points within the search radius are ranked by lowest to highest point and assigned a corresponding value between 0 and 100. This measure was intended to represent shrub shape and can particularly be useful for recognizing overhanging canopies. Also, characterizing shrub leaf and branch features, the angles between a point and all its neighbors within a search radius were extracted and averaged as the negative openness for each point. The negative openness refers to a conical view looking downwards, while oppositely the positive would be pointing up. It were originally developed for pixel‐based terrain modeling by taking the mean openness angle from the eight neighboring pixels (or more, depending on the search radius) in each cardinal and intercardinal direction (Yokoyama, Shirasawa, & Pike, 2002). However, all points in any possible direction within the search radius were considered for the 3D point attributes and therefore not necessarily restricted to eight directions.

From the initial set of 17 structural features, we retained seven after testing for autocorrelation and variable importance (Table 1). Variable importance was calculated manually based on the decrease in overall accuracy when leaving one variable out in the classification process (Appendix 2). Variables with high pairwise correlation (ρ > .75) and low importance scores were removed from the dataset.

TABLE 1

Variables selected for final classification

Variable	Overall importance score	Search radius (m)	Morphological features	Ecological relevance (presumed)	Computation description
Pseudowaveform	0.048	0.5	Shrub growth form	Increases with structural complexity (roughness)	Point height variances from a normal‐fitted plane
Point count	0.013	0.5	Shrub growth form	Increases with surface area and density	Number of points
Ground points count	0.006	0.1	Shrub growth form	Increases with light penetration and decreases with density	Number of points including ground points
Rank	0.002	0.1	Shrub shape and growth form	Highest values at upper surface and lowest at bottom canopies	Relative position in neighborhood. Max value = 100 and min value = 0
Negative openness	0.005	0.1	Shrub shape, leaf orientation, and branch arrangement	Increases with surface gaps	Angular distance between points from a nadir view
Point distance	0.025	0.1	Leaf structure and branch arrangement	Decreases with density	Average linear distance between points
No. of echoes	0.004	0.1	Leaf structure and branch arrangement	Decreases with light penetration	Average number of strongest and last return echoes

Machine learning for 3D classification of Cytisus scoparius and other shrub species

For the 3D classification of the shrub point cloud, we utilized the built‐in classification procedure in the OPALS software. It uses the tree‐based decision algorithm (De'ath & Fabricius, 2000) termed recursive partitioning via the rpart package for R (Therneau & Atkinson, 1997). The algorithm operates by dividing the dataset for each chosen variable separately, that is, it finds the splitting value for tree branches which results in the purest nodes, that is, most homogenous. A perfectly pure node refers to all observations being assigned the correct label caused by the split value. When the variables could not further increase the node purity, the resulting decision tree was pruned with the complexity parameter set to 0.001 for a final simplified decision tree. We fixed the complexity parameter after a trial process where we lowered the parameter stepwise until the accuracy started to decrease. The aim was to run the classification with the lowest possible complexity parameter.

The classification accuracy was assessed by randomly stratifying 90% of the reference data as training and 10% as testing data. The quality and reliability of accuracy assessments are affected by the reference data input and sampling strategy (Millard & Richardson, 2015). We therefore also assessed a 70% training to 30% validation data split (Appendix 2). Stratification was done among the defined classes. We generated 100 such 90/10 and 70/30 training and validation datasets as input for the OPALS classification algorithm. For each set, classification accuracy was assessed from the resulting confusion matrices using R to calculate overall accuracy and Kappa values, while the class‐wise accuracies were evaluated with precision, recall and the harmonic mean between the two, termed “F1” (Forman, 2003). The 70/30 split resulted in <0.5% decrease in overall accuracy and Kappa coefficient (Appendix 2). To increase the amount of training data used in the final model for extracting biomass metrics, we applied a model validated by the 90/10 split (Appendix 2). From the predictions, we obtained the class probabilities allowing fuzzy classifications, that is, the membership of a class is represented by a probability value between 0 and 1 rather than a Boolean value (true or false) as with traditional hard‐boundary classification(Foody, 1996; Zlinszky & Kania, 2016). Finally, a classification was performed with 100% of the reference data and printed into the full point cloud for visualization purposes.

Assessing variable and model transferability

We used a Wilcoxon signed‐rank test to evaluate the variable signatures for each shrub species in the validation samples from 2017 and 2018 (Appendix 2). Furthermore, the classification procedure was applied to a merged point cloud, consisting of the mixed signatures from 2017 and 2018.

Additionally, we applied a similar classification procedure as described above, but with the variables projected into 2D raster data (Appendix 2). The 2D classification was performed with R statistics 3.5.0 (R Core Team, 2016) allowing us to test alternative classification algorithms.

Biomass estimation and change detection

To estimate the relationship between actual measured biomass and structural information for each of the 10 harvested C. scoparius shrubs, we extracted digital volume metrics based on the NormZ variable equivalent to vegetation height (Appendix 3). The NormZ pixel values were computed using mean, max, and range (max–min) values from the points and thereby representing volume metrics by multiplying the height with the pixel area (25 cm2). We then used the manually constructed 3D point clouds from the field GNSS measurements to delimit and extract summed volume values for each of the harvested shrubs (Appendix 3).

The extracted volume metrics were all compared for correlation with nonparametric Spearman's rank correlation using R statistics 3.5.0 (R Core Team, 2016). Furthermore, we developed linear models from the volume metrics to explain biomass from the harvested samples and evaluated them by calculating adjusted R squared values ( ) via the lm() function in R. We used values from the models to cope with the relatively low sample size of 10 and thereby avoid making too optimistic conclusions. We applied the model coefficients (Figure 3) from the best fit (maximum NormZ) to calculate a biomass estimate for the rasterized variable across a 6.7 ha area and for both datasets (2017 + 2018). During the raster projection, we incorporated class probabilities from the final fuzzy classification models (see accuracy assessment in Appendix 2) to exclude points classified as C. scoparius with <60% probability for one set of maximum NormZ values. Likewise, a second set of values was extracted by adjusting this probability threshold more strictly to 80%. At last, the change in biomass was mapped and aggregated from 5 cm resolution to 2.5 m grid cells to emphasize the change.

FIGURE 3

Scatterplot with actual measurements of Cytisus scoparius biomass on Y‐axis and LiDAR‐based volume on X‐axis. Dashed blue lines represents corresponding linear regression models: Biomass ~ Volume for mean, max and range values. ( values for the models: mean, max and range values: .72, .77, and .61 with p‐values respectively: .0012, .0005, and .0048). Formulas: Y(Mean) = 92.49 + 4,575.88X; Y(Max) = 77.12 + 4,696.16X; Y(Range) = 92.98 + 12,327.87X

RESULTS

Detection of Cytisus scoparius in a 3D landscape of points

After georeferencing and noise filtering, our UAS LiDAR system generated 59 million and 53 million points during the 2017 and 2018 flight campaign, respectively, for the 6.7 ha area. For 2017, 12.6% of these points were classified into the shrub layer, and 6.1% in 2018. Likewise, 29.8% and 40.5% of the shrub layer points were classified as C. scoparius in 2017 and 2018, respectively (see Appendix 2). The overall accuracy from the 90/10 split classification was 86.9% for 2017 and 95.2% for 2018, while the Kappa coefficient resulted in 83.7% for 2017 and 92.9% for 2018. The 70/30 split resulted in <0.5% decrease in overall accuracy and Kappa coefficient (Appendix 2). Focusing exclusively on the C. scoparius class, the F1 measure was 96.2% in 2017 and 98.4 in 2018 (Table 2). The accuracy assessment from the specific model used for obtaining biomass differed with +0.2% in 2017 and +0.5% in 2018 from the mean of the iterative model validation with 90/10 splits (see details in Appendix 2).

TABLE 2

Accuracy assessment of the point cloud classifications of Cytisus scoparius from autumn 2017 and spring 2018. Classification accuracies are averages from 100 model iterations of randomly selected training/validation data (90/10% split) and with standard deviations. Overall accuracy and Kappa coefficient evaluate the classification of all classes, while the F1 score (Forman, 2003) assesses the performance in predicting the C. scoparius class. The last row presents the results from a merged classification including both the 2017 and 2018 point clouds

Time period	Overall accuracy		Kappa coefficient		F1 score—C. scoparius
	Mean	SD	Mean	SD	Mean	SD
Autumn 2017	86.9	0.006	83.7	0.007	96.2	0.005
Spring 2018	95.2	0.004	92.9	0.006	98.4	0.003
Merged	83.2	0.007	78.0	0.010	96.3	0.004

The analogous merged classification showed a decrease in overall accuracy (83.2%) and Kappa coefficient (78.0%), while the C. scoparius F1 remained equally high (96.3%) (Table 2). For the more traditional 2D classification approach, the random forest classifier performed in general better with overall accuracies of 45.6% for 2017 and 42.6% for 2018 (see Appendix 2). However, the recursive partitioning classifier did not perform markedly worse (overall accuracy of 47.1% and 38.9% for 2017 and 2018).

The color‐coded point cloud visualizing the resulting shrub classification (Figure 4a) revealed that C. scoparius shrubs were detected across the area in various densities, from single individuals to large thickets. When inspecting the classification probabilities in the point cloud (Figure 4b), we recognized a greater certainty to predict C. scoparius in the leaf‐off period in 2018 than in the leaf‐on period. Furthermore, it is noticeable that C. scoparius shrubs were detected underneath the forest canopy as well (Figure 4a1).

FIGURE 4

Classification results visualized directly in the 2018 point cloud in a vegetated area (a1) and in more open land (a2). The general classes are showing both, the main vegetation classes obtained from thresholding (Ground, Low vegetation, High vegetation) as well as the classes obtained by the tree‐based classification. The orange colors represent points classified as Cytisus scoparius. Lower figures show 3D clips of the 2017 (b1 to b3) and 2018 (b4 to b6) C. scoparius point clouds. Points are here color‐coded with the resulting class probabilities ranging from 60% (dark blue) to 100% (light green). The generally classified ground points are colored in grey

Cytisus scoparius biomass estimation

The LiDAR‐derived volume metrics of NormZ (average; maximum; range) correlated well with the biomass measurements of the 10 harvested C. scoparius shrubs (Spearman's ρ = .87; .88; .88). However, the inclusion of very small shrubs in the harvested samples challenged the LiDAR detection of biomass (Appendix 3). The maximum volume resulted in the best linear fit with  = .77 (Figure 3). The mean and range metrics performed worse in a linear model with  = .72 and  = .60, respectively. The accuracy assessment from the specific model used for obtaining biomass metrics is presented in Appendix 2.

Cytisus scoparius biomass change

We extracted the total biomass sum from the points with >60% and >80% probability of being C. scoparius. This resulted in 7,500.4 and 5,257.6 kg in 2017 and 2018, respectively, for >60% and, likewise, 5,320.9 and 4,993.0 kg for >80%, in an area of 6.7 ha. For the comparison, we included only the overlapping areas that had an average point distance below 3 cm. On the landscape scale, this resulted in an average biomass decline of C. scoparius from autumn 2017 to spring 2018 of 33.4 and 4.9 g/m2 for the 60% and 80% probability thresholds, respectively. However, on a local scale the distribution of biomass changes in the area varied but was similar for both thresholds. An upscaled visualization from 5 cm resolution to grid cells of 2.5 m × 2.5 m revealed a pattern of larger decreases in biomass to be identified in especially the northeastern part of the area, while in other parts, we observed no or a slight increase in biomass between the 2 years (Figure 5).

FIGURE 5

Orthophoto from the study area (Kortforsyningen Danmark, 2019) with overlay of change in maximum Cytisus scoparius volume between 2017 and 2018 UAS LiDAR survey. A C. scoparius class probability threshold is applied and only allowing points classified with >60% (top) and >80% (bottom) probability to be included for change detection. The change is upscaled from 5 cm to 2.5 m2 for visualization. Blue squares represent growth, while red squares show a biomass decrease. White squares indicate no or little change

DISCUSSION

With the sole use of UAS LiDAR‐derived structural information, we identified shrub classes with an overall mean accuracy of more than 86.9% regardless of the year and detected C. scoparius with at least 96.2% F1 accuracy in the point cloud (Table 2). Furthermore, using a simple volume metric (NormZ) from the classified C. scoparius point cloud, we explained 77% of the variation in actual harvested biomass (Figure 3) and quantified a reduction of 327.8 kg during the winter period from autumn 2017 to spring 2018, assuming the 80% probability threshold to be most accurate.

LiDAR‐derived structural features for shrub classification

Our study demonstrates that ecologically meaningful features can be extracted from UAS LiDAR point clouds which represent the structural appearance (i.e., growth form, shape, and leaf/branch orientation and arrangement) of different shrub species, especially C. scoparius. We targeted the variables to represent shrub species on varying levels of scale, ranging from small‐scale leaf characteristics to larger‐scale shrub growth form, density, and shape. For example, C. scoparius is characterized as being light competitive (Peterson & Prasad, 1998), hence, reducing the available light beneath the canopy. However, compared with, for example, the compact canopy of Juniperus communis, it may allow more light to travel through, especially during leaf‐off periods, eventually resulting in more ground hits of the LiDAR beams.

The variable that explained most of the variation, Pseudowaveform, was computed with a relatively large search radius (0.5 m) and therefore represents the variation in the structure of shrub branches and leaves (Table 1). By calculating on this relatively larger scale, we expected the variable to be less affected by differences in leaf abundance between the two flights. From analyzing the shrub signatures, we find Pseudowaveform to be stable between the two seasonal states and throughout the classes, except for Betula (Appendix 3). Structural complexity in vegetation is often assessed by the variance of derived LiDAR metrics relative to a horizontal plane (Kane et al., 2010; Kukunda et al., 2019). Textural metrics in image analysis can similarly be used to distinguish structurally dissimilar species (Oldeland, Naftal, & Strohbach, 2017), but are limited to detect only 2D surface differences. Here, we utilize information of 3D vegetation structure. The variance for example will be set relative to a normal‐based plane. Measurements of shrub growth form may be difficult to conduct and define in field‐based surveys; hence, there is a focus on retrieving metrics as, for example, diameter at breast height (DBH; Wieser et al., 2017). Nonetheless, the high importance of Pseudowaveform for classifying shrub species suggests that shrub growth form is highly relevant for distinguishing species. Ground points count and no. of echoes express the amount of light penetrating a shrub to represent density. Density is known to be highly variable among species (Van Leeuwen et al., 2011) while light penetration may also reflect habitat characteristics for the species composition (Moeslund et al., 2019). Density was also computed more directly with the amount of LiDAR points within a 50 cm search radius (point count), while point distance was computed on a smaller scale mainly emphasizing single leaf/branch structure and arrangement. Because of the direct influence from density, it remains challenging to estimate biomass in areas with high species heterogeneity (Wijesingha, Moeckel, Hensgen, & Wachendorf, 2019), and to improve accuracy, it is necessary to separate the species by a classification process. The negative openness was used as a metric for differentiating the structural complexity in especially the surface of shrubs. Specifically, the openness variable detects concaveness (Yokoyama et al., 2002) or in this case holes in the vegetation surface. From the shrub signatures, we see that negative openness is seasonally affected by a species like Betula (Appendix 3) showing lower concaveness in the leaf‐on than leaf‐off period. This can possibly be explained by leaves forming a relatively round and smooth canopy, whereas the branch gaps will be more exposed during leaf‐off period. Finally, to differentiate between C. scoparius and low tree species in the shrub point cloud with overhanging canopies, for example, Malus sp. or Prunus cerasifera, the classification benefitted from the Rank variable. Here, lower points from a hanging branch without points underneath would be ranked lower than C. scoparius branches, which often will have lower points beneath, from low vegetation or ground. Thus, computing ecologically meaningful variables from 3D point clouds is possible but needs to be based on ecological knowledge and targeted on the focal species or vegetation class.

The performance of the LiDAR‐derived structural variables to detect and distinguish shrub species in this study demonstrates a promising use at the level of detail obtainable with a drone‐based platform (Table 2). In addition to the increased spatial resolution and structural detail, our landscape study of C. scoparius differs from Hill et al. (2016) in being independent from observations during flowering periods. However, phenological events can also be considered an important aspect of understanding the development of shrub encroachment. In particular, an increase in temporal resolution would be beneficial for remote sensing studies of vegetation to understand how such variables vary with seasonality. A study by Müllerová et al. (2017) demonstrates an example of this, by recognizing two invasive species from differently scaled images and throughout the season. Hence, using spectral information in combination with LiDAR could extend the monitoring possibilities even further.

Workflow and classification challenges

The use of LiDAR‐derived features poses many challenges and demands a novel way of developing and understanding these 3D measures in an ecological context. Thus, further improvements are possible by implementing ecological or biological knowledge in the computed variables. Depending on the focal species, the LiDAR‐derived features can be adjusted to fit specific morphological characteristics, and a more general approach could be developed to assure transferability for many shrub species. For C. scoparius, the change in structural signature between leave‐on and leave‐off is expected to be low due to the small leaves (Appendix 2). However, with the onset of flowering, this might change for some of the variables. Alternatively, for other shrub species it might be important to develop variables that are independent of seasonal states to assure transferability of the classifier. For C. scoparius, we found a significant change in the sample distributions in three of the density variables (PCount_ground, PCount, and PDist) which indicates that seasonal variability might be seen in these variables. However, it is the combined signature that is important for the classification process and the distinction among species. When evaluating normalized vegetation height (NormZ) used for biomass estimation, we find no difference between the two flight dates in C. scoparius, which is different for Betula species. These findings highlight the need for variable adjustments and/or seasonal timing, to obtain reliable biomass estimations for a given target species.

We implemented a workflow to improve the assessment of shrub biomass during encroachment. During the development, we experienced several challenges which could affect the results and need careful consideration and evaluation during implementation: (a) misalignment reduce spatial accuracy, (b) thresholding of major vegetation classes, for example, trees versus low vegetation, (c) sampling of training data, and (d) selection and development of targeted structural features. Below we outline these challenges and provide perspectives on future directions.

Because of observed misalignments in the combined point clouds, resulting from different laser beams, flight strips, and uncertainties in absolute GNSS positioning, we rectified the relative accuracy in the datasets, which negatively influence the absolute accuracy, that is, absolute geographic position of LiDAR points. However, the decrease in precision of XYZ‐coordinates was relatively small (RMSE: vertical <1.4 cm and horizontal <5.1 cm). When considering the scale of which we recognize C. scoparius, we believe that a high relative accuracy among points is more important for describing the structural features of a shrub properly. However, in monitoring programs, fixed features (i.e., ground control points) allowing to align multiyear data would be advisable.

Threshold approaches to filter data noise as well as for doing a first general classification are essential for effective data processing. However, finding meaningful thresholds in heterogeneous plant communities with highly variable structures is difficult and hard to generalize across a larger area. Therefore, deciding on exact thresholding values often is a trade‐off issue, by either including an amount of noise or excluding useful data. In our study, we emphasized to include as much data as possible, accepting some levels of uncertainty. Still, with broadly defined threshold values, points representing a target species may disappear in these first steps of filtering and classifying. For example, we attempted to include training data from the small‐sized shrub Calluna vulgaris; however, these points were set aside in the general class “Low vegetation” and thus not included in the tree‐based classification. Reducing the dataset is crucial for optimal processing speed and point cloud visualization but choosing general classes and thresholds must be evaluated and depends on the target species. If targeted, the C. vulgaris shrubs could most likely be identified by running the classification procedure on the low vegetation points and could perhaps be implemented in a future study including herbaceous vegetation.

Generally, one should aim for balanced training/validation data (Millard & Richardson, 2015). However, this can be challenging when working in natural areas with a patchy and heterogeneous distribution of vegetation. First, a random sampling strategy to collect data objectively is desired, a strategy we pursued with the transect approach. However, this method may only cover the more dominant species and is not enough to catch the full variation of shrub species in an area. Second, it might be impossible to equally sample rare species occurring only in few places, even if all individuals are sampled. We therefore included datapoints outside our transects from less abundant shrub species, as it would not present the full picture of the area if ignoring those. When covering extents of several hectares, we are aware that some species are possibly missing in the training data, but that does not mean that they do not exist in the point cloud.

It is important to design point‐based attributes to emphasize structural features, corresponding to the species or vegetation type of interest. The described workflow is semiautomatic, given a comprehensive input dataset of UAS LiDAR point clouds, training data for classification, and ground truth data for developing the biomass model. Depending on the study focus and target species, it is possible to design point‐based attributes that emphasize structural features, corresponding to a species or vegetation type of interest. For example, different search radii are required when studying differences in trees or herbs matching the level of scale. Likewise, attributes can be modified to detect specific structural features of species or individuals, for example, deformity caused by pathogens. The structural uniqueness, as defined by growth form, density, roughness, and light penetration seen in different shrub species may be less influenced by seasonality and physiological stress factors (i.e., water and nutrient availability) than spectral signatures. However, when attempting to distinguish structural similar species including spectral information in the classification procedure might help (Zlinszky et al., 2014).

Spatial patterns of C. scoparius biomass change

Depending on the probability threshold for detecting a C. scoparius shrub, we estimated a decrease in biomass change of 4.9–33.4 g/m2 between the two mapping dates. Nevertheless, we found a consistent spatial pattern of biomass change for both thresholds, indicating a large decline in the northeastern part of the area and no or a slight increase in the southwestern part (Figure 5). There might be several reasons for this observed spatial pattern:

First, the slight increase in biomass could result from late autumn or early spring growth. The shrubs are not expected to grow during a Danish winter season; however, the timing of flight campaigns may still include some late autumn or early spring growth. In addition, C. scoparius branches contain chlorophyll, which might be able to induce growth in warmer periods and favor early growth in spring. The reason why biomass increase was mainly observed toward the southwest of the study area could be low inter‐ and intraspecific competition or variation in microclimatic conditions. A freestanding individual without neighbors is more likely to grow and expand due to lack of competition when the conditions are suitable. In contrast, competing neighbors may supress growth of C. scoparius individuals in densely populated surroundings. The possible intraspecific interaction has previously been raised in the literature (Paynter, Fowler, Memmott, & Sheppard, 1998), and in accordance with this, the observed spatial pattern in biomass changes suggests that C. scoparius in the open land (midwestern part) tended to slightly increase, while the denser stands toward the northeast of the study area showed a decline in biomass (Figure 5).

Second, the reasons for a biomass decline could be either leaf fall (Peterson & Prasad, 1998), which due to the summed 2.5 m grids, is expected to be larger in denser stands, or external factors, such as grazing or frost. Galloway cattle and Exmoor ponies graze the area as part of a rewilding initiative, and no additional food sources are supplied to the animals. Hence, at wintertime when the green vegetation becomes sparse, the animals might feed on hardy shrubs like C. scoparius, or cause damage to shrubs by trampling. Therefore, one alternative hypothesis explaining the spatial pattern of C. scoparius biomass change is that the animals are favoring the northeastern part of the study area. This area is substantially more forested, and the grazers might benefit from shelter provided by the trees during harsh winter conditions. Alternatively, abiotic conditions linked to topographic features, such as light availability and freezing temperatures during winter might limit plant growth and harsh winters would potentially cause C. scoparius to die (Peterson & Prasad, 1998). Extracting such topographic variables describing light availability and protection against freezing could provide predictors for the observed change in biomass. The availability of terrain topography data from the UAS LiDAR point cloud would facilitate such a follow‐up study at very high resolution.

While the observed pattern in biomass change is independent of the probability thresholds for C. scoparius detection, the absolute biomass change might be overestimated. When comparing C. scoparius class probabilities in three areas with varying population densities (Figure 4b), there are more C. scoparius points in 2017 than in 2018, especially in areas with dense stands. In the leaf‐on period (2017), a larger number of forbs, tall grasses, and shrubs have probably been included in the training data for C. scoparius, as they grow within and around shrubs. This probably also explains why in 2017 more C. scoparius points have been classified with lower probabilities than in 2018 (Figure 4b and Appendix 2). More LiDAR points included from entangled vegetation can lead to a larger volume estimate per individual and probably cause slight overestimation of C. scoparius biomass in 2017.

Relevance for nature management

The spatial variation in C. scoparius biomass change indicates that shrub dynamics differ in the subareas of the mapped area, possibly due to varying importance of ecological drivers, also yet to be studied (Figure 5). Our findings highlight that detection of change in shrub density and biomass with high resolution is important when assessing shrub encroachment in monitoring programs for nature management. Traditionally, remote sensing studies have classified vegetation with a two‐dimensional approach from spectral information or rasterized LiDAR information. These raster‐based methods are efficient in identifying cover or presence of certain grassland species and vegetation types (Hellesen & Matikainen, 2013; Zlinszky et al., 2014). However, to quantify biomass changes and to fully understand the effect of shrub encroachment on plant diversity more comprehensive knowledge is needed. While LiDAR observations from a manned aircraft were applied to map coarse‐scaled shrub encroachment characteristics from a single species (Sankey et al., 2013), we are now able to separate species or genera directly in the point cloud and to detect fine‐scale biomass dynamics from a target species by utilizing a drone platform. With the use of the established workflow based on a UAS LiDAR system, we provide a new approach for monitoring shrub species dynamics. In this study, C. scoparius is covering the spatial extent corresponding to management operations, but also capturing the local‐scale information needed for detecting change and its spatial variation. Again, the benefits of using LiDAR are the vegetation penetration ability (Lefsky et al., 2002) and our findings suggest that LiDAR‐derived point clouds are of such a quality, that detecting species of interest even beneath a covering canopy is achievable, as, for example, forest understory species. This will make it possible to monitor and guide management programs for noxious invaders such as Rhododendron ponticum (Sanders, 2017) or rare species such as Allium tricoccum beneficial for indicating favored nature or ecological conditions (Leduc & Knudby, 2018). In a study by Chance et al. (2016), LiDAR is utilized to map the distribution of two invasive shrub species in an urban environment with reported lower accuracies for shrubs beneath closed canopies. The importance of mapping understory vegetation is further highlighted in a review by Hernandez‐Santin, Rudge, Bartolo, & Erskine, 2019, and the findings from our study should be relevant and encouraging to further studies in this direction.

Shrub encroachment is influenced by several factors that makes it a complex issue for management and conservation at landscapes and regional scales (Cao et al., 2018). The presented workflow provides a potential improvement on several aspects of nature monitoring. By the semiautomatic processing and the use of a UAS‐based platform, the temporal resolution is no longer limited to snapshot observations, for example, once a year or less. Hence, field operations can be carried out on multiple times through a season, with the inclusion of beneficial seasonal states of vegetation. However, we also demonstrate that UAS LiDAR data allow for comparable measures of shrub encroachment even in different seasonal states and hereby displaying the potential multipurpose use from this type of data. Our comparison with a similar 2D remote sensing classification approach (Appendix 2) highlights the benefits from using three‐dimensional information. Hence, we increase the spatial dimensionality from traditional 2D mapping studies by recognizing shrub taxa in three‐dimensional space contributing to reliable biomass estimation which is much needed for deeper understanding of shrub encroachment dynamics. With the fine‐scale information about biomass, it will be easier to detect and manage dense stands of C. scoparius or other shrub species during encroachment with potential negative effects on plant diversity (Kesting et al., 2015).

CONCLUSION

Our study presents a novel method for assessing shrub dynamics in 3D, based on the arrangement and orientation of points in space. We demonstrate an efficient way of determining specific structural features for classifying shrub species. Seven different point‐based structural variables were developed based on ecological knowledge to distinguish C. scoparius from other vegetation. Using a tree‐based classification procedure, we identified 11 different shrub species with an overall accuracy of 86.9% and 95.2% in two independent point clouds acquired in 2017 and 2018, respectively. Derived from the point cloud derived height‐based maximum volume of C. scoparius shrubs, we established a linear model to explain 77.4% of the variation seen in actual measurements biomass from 10 harvested and weighed individuals. Projected biomass in C. scoparius with an 80% class probability declined by 4.8 g/m2 between 2017 and 2018. However, we found substantial spatial variation in biomass change with a patchy distribution of C. scoparius decrease and or growth. These findings suggest that the combined information from small‐ to large‐scale observations is necessary to fully understand shrub encroachment and highlights the potential application of this method for fine‐scale management of C. scoparius populations on a landscape scale. Moreover, our study encourages more studies investigating the causes of C. scoparius biomass change variation and its effect on plant diversity.

Our study contributes with new and advancing methods within the field of remote sensing, ecology and nature management by demonstrating an application for ecological monitoring. Natural systems in Denmark, including grasslands, are monitored by the Danish national monitoring program (Svendsen, Bijl, Boutrup, & Norup, 2005). The promising results of detecting change in biomass related to encroachment of C. scoparius may inspire improving some of the traditional field‐based measurements toward a more objective evaluation on the state of natural systems. The presented work potentially improves monitoring programs in nature management as it allows to quantify biomass change timely and at a fine spatial scale informing managers to implement management strategies that help sustaining biodiversity during shrub encroachment.

CONFLICTS OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Bjarke Madsen: Conceptualization (lead); data curation (lead); formal analysis (lead); investigation (lead); methodology (lead); project administration (lead); resources (lead); visualization (lead); writing – original draft (lead); writing – review and editing (lead). Urs A. Treier: Conceptualization (lead); data curation (equal); formal analysis (equal); investigation (equal); methodology (equal); project administration (equal); resources (equal); supervision (equal); validation (supporting); visualization (supporting); writing – original draft (supporting); writing – review and editing (supporting). András Zlinszky: Formal analysis (supporting); investigation (supporting); methodology (supporting); software (supporting); validation (supporting); writing – original draft (supporting); writing – review and editing (supporting). Arko Lucieer: Investigation (supporting); supervision (supporting); validation (supporting); writing – original draft (supporting); writing – review and editing (supporting). Signe Normand: Conceptualization (equal); formal analysis (equal); funding acquisition (lead); investigation (equal); methodology (equal); project administration (equal); supervision (lead); validation (equal); writing – original draft (supporting); writing‐review and editing (supporting).

TABLE A1.1

Calculated root mean square error (RMSE) and mean absolute difference (MAD) from vertical deviation measurements in CloudCompare before and after relative adjustment (*_adj)

	Vertical deviation				Horizontal deviation
	2017	2018	2017_adj	2018_adj	2018	2018_adj
RMSE	0.039484	0.060631	0.076016	0.091212	0.014	0.064713
MAD	0.018583	0.025667	0.03175	0.030532	0.00696	0.03888

TABLE A2.1

Variable importance based on overall accuracy, completeness and correctness scores. The three measures were summed without additional weights (sum) and ranked from best to worst (1–17) according to each measure before averaging (rank). Color indicates the variables that are kept in the models for the next step (green), and those left out due to lack of importance with a rank threshold ≤ 10 (Red) or autocorrelation with pairwise cor >.75 (yellow)

Variable	Overall accuracy	Completeness	Correctness	Rank	Sum
PseudowaveL	0.047808195	0.00902837	0.047616142	1	0.104453
PDist	0.025160961	0.00811807	0.000653427	5.333333	0.033932
PCountL	0.013440036	0.002558371	0.01761165	3.666667	0.03361
Pcount_ground	0.006272326	0.002858146	0.008466972	5	0.017597
NoOfEchos	0.003692739	0.004184393	0.006445916	6	0.014323
NegOpen	0.004667219	2.81416E−05	0.009419098	7	0.014114
Range	0.004524673	−0.005458463	0.014149932	9	0.013216
NormZ	0.006381884	0.003413731	−9.15274E−05	7.333333	0.009704
Planarity	0.001868897	−0.000228497	0.005266631	10	0.006907
PCountM	0.003716534	−0.000250142	0.003000489	10	0.006467
Rank	0.001781419	0.002352647	0.001707153	10	0.005841
PseudowaveS	0.00252669	0.000889567	0.001788535	9.333333	0.005205
Linearity	0.000436458	0.000820137	0.000593598	12.33333	0.00185
PDensL	0.001645286	0.000281219	−0.000591189	12.66667	0.001335
PCountS	0.001047605	−0.002076538	−0.001163251	14.66667	−0.00219
StdDev	−0.001450466	−0.005435544	0.004469112	13.33333	−0.00242
PosOpen	−0.003849571	−0.003340636	−0.004994585	16.33333	−0.01218

TABLE A2.2

Differences in shrub signatures between 2018 and 2017 flights for each variable and corresponding Wilcoxon p‐values comparing the samples. p‐Value significance codes are illustrated by circles: p‐val ≤.05 (black); p‐val ≤ .1 (grey); p‐val >.1 (white)

TABLE A2.3

Accuracy measures from the final models used for calculating biomass from the 2017 and 2018 flights. Accuracies derived from the resulting confusion matrix are presented as overall accuracy and Kappa coefficient. The class‐wise accuracies are presented as precision, recall, and F1 measures

Accuracy assessment—Final model
	2017	2018
Overall Accuracy	0.8706	0.956609
Kappa	0.8394	0.935633

TABLE A2.4

Accuracy measures from the merged classification including the 2017 and 2018 point clouds. Accuracies derived from the resulting confusion matrix are presented as overall accuracy and Kappa coefficient. The class‐wise accuracies are presented as precision, recall, and F1 measures

Accuracy assessment—Final model
	70 split	90 split
Overall Accuracy	0.829	0.832
Kappa	0.775	0.780

TABLE A2.5

Overall accuracy and Kappa values resulting from a 2D raster classification of shrub species. Both datasets from 2017 and 2018 were classified with four different classifiers (Model name)

Model name	Overall accuracy	Kappa
Rpart17	0.4706	0.1915
Rpart18	0.3889	0.1449
RF17	0.4559	0.2465
RF18	0.4259	0.1838
GBM17	0.3971	0.1761
GBM18	0.3889	0.1518
DNN17	0.3971	0
DNN18	0.4815	0

TABLE A3.1

Values from the 10 harvested individuals by measured dry weight biomass in gram (Weight) and extracted volume in m3 from exported raster data (Mean, Max and Range)

Individual	Weight	Mean	Max	Range
1	132	0.01196	0.00786	0.00033
2	430.8	0.07370	0.09342	0.03822
3	626.5	0.13262	0.11979	0.02830
4	33.7	0.00550	0.00555	0.00010
5	104.6	0.01923	0.01183	0.00256
6	287.6	0.03568	0.04844	0.02597
7	48.8	0.00131	0.00131	0.00000
8	153.1	0.01802	0.02552	0.01223
9	554.6	0.03844	0.04625	0.01692
10	107.1	0.00313	0.00365	0.00102