Intra-Annual Variabilities of Rubus caesius L. Discrimination on Hyperspectral and LiDAR Data

The study was focused on a plant native to Poland, the European dewberry Rubus caesius L., which is a species with the ability to become excessively abundant within its original range, potentially causing significant changes in ecosystems, including biodiversity loss. Monitoring plant distributions over large areas requires mapping that is fast, reliable, and repeatable. For Rubus, different types of data were successfully used for classification, but most of the studies used data with a very high spectral resolution. The aim of this study was to indicate, using hyperspectral and Light Detection and Ranging (LiDAR) data, the main functional trait crucial for R. caesius differentiation from non-Rubus. This analysis was carried out with consideration of the seasonal variability and different percentages of R. caesius in the vegetation patches. The analysis was based on hyperspectral HySpex images and Airborne Laser Scanning (ALS) products. Data were acquired during three campaigns: early summer, summer, and autumn. Differentiation based on Linear Discriminate Analysis (LDA) and Non-Parametric Multivariate Analysis of Variance (NPMANOVA) analysis was successful for each of the analysed campaigns using optical data, but the ALS data were less useful for identification. The analysis indicated that selected spectral ranges (VIS, red-edge, and parts of the NIR and possibly SWIR ranges) can be useful for differentiating R. caesius from non-Rubus. The most useful indices were ARI1, CRI1, ARVI, GDVI, CAI, NDNI, and MRESR. The obtained results indicate that it is possible to classify R. caesius using images with lower spectral resolution than hyperspectral data.


Introduction
Growing global phenomena such as land-use changes or habitat fragmentation, and the accompanying climate change yielding changes in the ecological and geographical ranges of species, lead to biodiversity loss [1][2][3][4][5]. Plant species that spread on a massive scale beyond their original geographical ranges and native species becoming excessively abundant within their original ranges often cause significant changes in ecosystems, including biodiversity loss [6][7][8].
The present study is focused on a Rubus caesius L. species from the brambles genus (Rubus), which has traits that favour spreading, such as low trophic requirements and quick adaptation to changing habitat conditions. In addition, by competing with other resident species, brambles limit the populations of other species by competition and strongly modu-strongly modulate species interactions and community composition, leading to a decline in biodiversity in the affected habitats [9][10][11][12], including protected in Europe Natura 2000 habitats.
Rubus caesius L. (European dewberry, R. caesius) belongs to the large and diverse genus of Rubus (bramble) within the Rosaceae family, covering about 430 and 750 species [13]. Numerous species and infrageneric taxa have been subsequently recognised, recently reaching about 1500 species worldwide [14], depending on the taxonomic approach. This plant is a half-bush reaching a height of 0.5-2 m with shoots lying on the ground. The main identification features of this species include the morphology of its stems (usually slender and pruinose, with short, slender, and needle-like prickles), leaves (three-foliate with broad-based, gibbous lateral leaflets; the stipules are broadly lanceolate) and generative organs: flowers (white flowers, appearing both on previous-year and current-year twigs, gathered in corymbose inflorescences) and fruit (ripe black drupeles pruinose) [15].
Rubus caesius grows in forest and shrub communities with natural features (mostly in riverside riparian woodland), as well as at the edges of forests but much more often in meadows (formed in the habitats of ancient riparian forests) and ruderal communities ( Figure 1) [15][16][17]. Dewberry is also a threat to non-forest natural habitats [16]. In many cases, the plant creates large, homogeneous patches. It also may cause problems as a weed, e.g., in maize and winter wheat [18,19] or other cereal crops (own observations). R. caesius depends mainly on vegetative reproduction using the buds on its roots and creeping stems and less on its seeds. Its developmental phases correspond to the seasons of the year, but R. caesius tends to flower throughout much of the growing season [20]. Thus, vegetative growth occurs from early spring, with the flowering and fruiting phase ranging from May to September. The first-year canes or shoots of R. caesius are pedunculated, round, straight, thin, glabrous, or slightly pubescent, waxy frosted, herbaceous bluish, glaucous, and tinged red where exposed to sun. Two-year-old shoots are partly woody at the base, branched, bloom and bear fruit, and usually die after two seasons. The leaves contain tannins, flavonoids, anthocyanin dyes, organic acids, vitamin C, vitamin B, vitamin P, inositol, and pectin [21]; moreover, Rubus caesius leaf extracts have medicinal value [22,23]. A detailed description of species variability during the growing season is described in Section 3.2.
Rubus caesius is widespread throughout Europe and Western and Central Asia [15] but is also observed outside its original range in North America [24]. In methodological guides prepared for the purposes of monitoring of specific protected ecosystems, Rubus is listed among the expansive species that threaten the following non-forest Natura 2000 R. caesius depends mainly on vegetative reproduction using the buds on its roots and creeping stems and less on its seeds. Its developmental phases correspond to the seasons of the year, but R. caesius tends to flower throughout much of the growing season [20]. Thus, vegetative growth occurs from early spring, with the flowering and fruiting phase ranging from May to September. The first-year canes or shoots of R. caesius are pedunculated, round, straight, thin, glabrous, or slightly pubescent, waxy frosted, herbaceous bluish, glaucous, and tinged red where exposed to sun. Two-year-old shoots are partly woody at the base, branched, bloom and bear fruit, and usually die after two seasons. The leaves contain tannins, flavonoids, anthocyanin dyes, organic acids, vitamin C, vitamin B, vitamin P, inositol, and pectin [21]; moreover, Rubus caesius leaf extracts have medicinal value [22,23]. A detailed description of species variability during the growing season is described in Section 3.2.
Rubus caesius is widespread throughout Europe and Western and Central Asia [15] but is also observed outside its original range in North America [24]. In methodological guides prepared for the purposes of monitoring of specific protected ecosystems, Rubus is listed among the expansive species that threaten the following non-forest Natura 2000 habitats: xerothermic grasslands Festuco-Brometea (code 6210) [25], mountain yellow Trisetum and Bent-grass hay meadows (Polygono Trisetion and Arrhenatherion) (code 6520), lowland hay Remote Sens. 2021, 13, 107 3 of 22 meadows (Arrhenatherion) (code 6510) [26,27], and also species-rich Nardus grasslands on siliceous substrates in mountain areas (and submountain areas, in Continental Europe) (code 6230) [27] and dry heath communities Calluno-Genistion, Pohlio-Callunion, Calluno-Arctostaphylion (code 4030) [28]. In the case of the meadows, the presence of Rubus most often testified to the absence of mowing. In the Natura 2000 European Ecological Network, invasive alien species (IAS) is one of the key parameters that provides information about the level of status conservation of every habitat [29]. Hence, quick identification of the presence of this species is essential for the effective protection of many natural habitats. To support assessments of the effects of global change on ecosystems, an important task is to develop easy methods and tools for data integration, processing, storage, and analysis. Detecting the species is the first step in monitoring and should be done as simply as possible.
For the Rubus genus, different types of data were successfully used for classification. Bramble R. fruticosus was effectively monitored using HyMap hyperspectral data with high accuracy (the best Overall Accuracy (OA) was 92%; kappa 0.715) [30]. Data were acquired in November (spring in Australia). HySpex images were used to classify three different expansive species including Rubus spp.; the F1 for the best dataset was equal to 0.97 for the Support Vector Machine classifier and 0.95 for Random Forest [31]. The Himalayan Giant bramble Rubus armeniacus was successfully classified (OA from 61.8% to 76.4%) using spectral bands from the Compact Airborne Spectrographic Imager (CASI) and Light Detection and Ranging (LiDAR) [32]. The same spectral CASI data combined with LiDAR-derived layers were used in the classification, where the OA was equal to 87.8% and the kappa coefficient was 0.75 [33]. The species Rubus caesius was successfully classified on the basis of the CASI and LiDAR data, where the User Accuracy (UA) for the class was 82%, and the Producer Accuracy (PA) was 58% [34].
Some of the multispectral satellite data were also tested for Rubus mapping. Satellite Landsat 8 and Sentinel-2 were used to identify the American bramble (Rubus cuneifolius) [35]. As a result, Sentinel-2 out-performed Landsat 8 in all seasons, with summer showing the highest accuracy (OA 77%, PA 80%, UA 45%). The Sentinel-2 bands of NIR, red-edge, and SWIR were crucial for increasing mapping accuracy. The images from SPOT 5 acquired in different seasons were used to identify Rubus cuneifolius as an invasive species in northern KwaZulu-Natal [36]. Different classifiers and approaches were tested, and the best OA for Rubus was 82% with a kappa of 0.71.
Based on this information, dewberry can be successfully identified, but it is necessary to collect optimal data. The high spectral and radiometric resolution makes it possible to identify brambles with high accuracy, but it is also possible to map large patches of Rubus using satellite Sentinel-2 images. For high spectral resolution data, much information must be processed over a long time using dedicated algorithms. Identification of the species through monitoring must also be precise and as fast and easy as possible. Thus, it is crucial to define the remote sensing data for proper identification.

Aim of the Study
Monitoring vegetation areas requires species identification that is fast, reliable, and repeatable. Dewberry classification is possible, which was proven in the studies mentioned above. At the same time, it has not yet been determined which remote sensing data are necessary for proper classification. The aim of this study is to define the following: • What type of remote sensing data differentiates R. caesius from non-Rubus depending on its coverage and the date used for image acquisition; • Which remote sensing data (spectral bands, calculated indices, and structural metrics derived from ALS) are the most discriminant of R. caesius under different growth and pigmentation phases.
To find the discriminating layers within different datasets, hyperspectral images and products of Airborne Laser Scanning (ALS) statistical analysis were used along with Linear Discriminate Analysis (LDA) and different ANOVA tests. These previously successfully used methods to discriminate vegetation class data and identify the most differentiated Remote Sens. 2021, 13, 107 4 of 22 layers within the spectral bands [37,38] or within the vegetation indices [39] for crop identification. Analyses were also used for analysis of different high-Arctic plants [40].

Airborne Data Acquisition
The study was based on hyperspectral and ALS data. Three instruments, two HySpex scanners (VNIR-1800 and SWIR-384), and an ALS scanner (Riegl LiteMapper) were placed in the same aircraft and combined as a HabitARS platform [41]. The data were acquired simultaneously: hyperspectral images in 451 spectral bands were acquired using a VNIR-1800 (range 400-1000 nm, 182 bands, spectral sampling 3.26 nm, spatial resolution 0.5 m), a SWIR-384 (range 950-2500 nm, 288 bands, spectral sampling 5.45 nm, spatial resolution 1 m), and ALS data acquired with a density of 7 points/m 2 [41]. Overflights were made during three campaigns in 2016: early summer (June 22)-C1, summer (July 24)-C2, and autumn (September 11)-C3. In each of the campaigns, data for the same area were obtained.
A total area of 307 ha was selected for the study, located in the Soła river valley near Oświęcim (southern Poland), where R. caesius occurs with high frequency and density ( Figure 2). The area was abandoned by agriculture (no treatments, including mowing).

Acquisition and Preprocessing of Ground Reference Measurements
The field study was done over three campaigns corresponding to the overflights: To establish R. caesius polygons, during the first campaign, 50 patches in which this species dominated with different percentages of coverage acquired every 10% (from 40 to 100%) were located, and within those patches, reference polygons were established ( Figure 3). Each polygon was located within a compact and homogeneous dewberry patch. Species coverage was assessed using the modified Braun-Blanquet method [42] in increments of 10%, where, e.g., 10% means that the coverage of this species was estimated between 5 and 15%. This method is based on field estimation of the percentage coverage of aboveground shoots in the patch and is commonly used in vegetation studies. The information about R. caesius coverage in reference was acquired in each campaign.
s. 2021, 13, x FOR PEER REVIEW 6 of 23 produced flowers or fruits; leaves and shoots herbaceous bluish, glaucous and tinged red where exposed to the sun glandulifera, Artemisia, Urtica) exceed the height of the R. caesius patches (1-1.5 m); plants of some species bloom profusely (Impatiens glandulifera)

C3
Mostly fruiting phase; fruits not massproduced, did not stand out from the dominant leaves; leaves and steams herbaceous bluish, glaucous and tinged red where exposed to the sun (more intense process than in C2) Discolouration of leaves of many species (both herbaceous and woody) appear; visible senescence of spring and early summer species plants (drying out) All polygons were processed before further analysis. For each campaign, we chose polygons without shadows that had not been mowed or grazed. For each campaign, based on visual interpretation, we also drew 30 polygons for trees.

Airborne Data Processing
The ALS point cloud orientation was carried out using RiProcess software [43]. The accuracy of this process was assessed at the level of 1 Σsigma = 0.010 cm; then, the data were automatically classified and corrected manually in the TerraSolid software [44]. The classified points were used to generate the Digital Elevation Model (DEM) and Digital Surface Model (DSM) with a 1 m resolution.
The first step in the processing of hyperspectral data was radiometric calibration (performed in the HySpex RAD program, Norsk Elektro Optikk AS, Skedsmokorset, Norway [45]) followed by parametric geometric correction of each pixel (based on ALS' DSM and Inertial Navigation System data) [41]. This process was performed in the PARGE en- Rubus caesius patches varied from single clumps to patches covering several square meters to several dozen square meters. No dry shoots or dry biomass were visible from above in any campaign. An observation chart was created for each of the polygons, which was completed during the second and third campaigns. A summary of R. caesius's developmental phases is included in Table 1.
During the three campaigns, we acquired a total of fifty polygons with vegetation different than R. caesius. These polygons covered the entire variability of the surveyed area, including plant communities in which R. caesius was recorded or dominated by species with similar morphologies. The following species and types of vegetation were present within the non-Rubus reference plots: Artemisietea vulgaris, Molinio-Arrhenatheretea, Phragmitetea and Cirsium arvense, and Solidago or Tanacetum. The affiliation polygons to the non-Rubus subclasses did not change during the vegetation season, so the class obtained once was constant throughout the year (from C1 to C3). Flowering and fruiting phase; both phases extended over time, not mass-produced flowers or fruits; leaves and shoots herbaceous bluish, glaucous and tinged red where exposed to the sun Some of the herbaceous plants (Solidago, Phragmites, Tanacetum, Impatiens glandulifera, Artemisia, Urtica) exceed the height of the R. caesius patches (1-1.5 m); plants of some species bloom profusely (Impatiens glandulifera)

C3
Mostly fruiting phase; fruits not mass-produced, did not stand out from the dominant leaves; leaves and steams herbaceous bluish, glaucous and tinged red where exposed to the sun (more intense process than in C2) Discolouration of leaves of many species (both herbaceous and woody) appear; visible senescence of spring and early summer species plants (drying out) All polygons were processed before further analysis. For each campaign, we chose polygons without shadows that had not been mowed or grazed. For each campaign, based on visual interpretation, we also drew 30 polygons for trees.

Airborne Data Processing
The ALS point cloud orientation was carried out using RiProcess software [43]. The accuracy of this process was assessed at the level of 1 Σsigma = 0.010 cm; then, the data were automatically classified and corrected manually in the TerraSolid software [44]. The classified points were used to generate the Digital Elevation Model (DEM) and Digital Surface Model (DSM) with a 1 m resolution.
The first step in the processing of hyperspectral data was radiometric calibration (performed in the HySpex RAD program, Norsk Elektro Optikk AS, Skedsmokorset, Norway [45]) followed by parametric geometric correction of each pixel (based on ALS' DSM and Inertial Navigation System data) [41]. This process was performed in the PARGE environment (ReSe-Remote Sensing Applications, Wil, Switzerland) [46]. Orthorectification was done using DSM from ALS data. As a result, the pixel location accuracy of RMS = 0.77 m was obtained as compared to orthophotomap with 0.1 m GSD. The atmospheric correction was done in the ATCOR4 program using the MODTRAN5 model (ReSe Applications GmbH, Wil, Switzerland) [47]. The final stage was a mosaic of series made in OrthoVista [48]. Next was the elimination of bands, where the absorption of radiation by water occurs, i.e., 21 bands were removed from the scope of the final scanner due to the range of absorption of radiation by water vapour in the atmosphere (above 2396.44 nm). Finally, an image consisting of 430 bands in the range of 416.18-2396.44 nm was obtained.
To assure the geometrical co-registration, HySpex mosaic and raster derivative products from ALS were integrated into a single set with a standard coordinate system, matrix start point, and pixel size equal to 1 m.
Hyperspectral data were used to calculate 49 remote sensing vegetation indices (VI) in the ENVI 5.3 program using the Spectral Indices function. The list of the indices with the equation is in Table S1.
Based on ALS data, raster layers were created from a point cloud using the Boise Center Aerospace Laboratory LiDAR (BCAL) tool [49]. Then, products from the Intensity and Vegetation groups were calculated. The first group of products obtained refers to the intensity of reflectance, while the second refers to the parameters related to the height of the plant cover. The Vegetation Products were calculated based on elevation and the maximum height of canopy (first return of the laser pulses). In addition, from the point cloud, products related to the vegetation structure were generated in the Orientation and Processing of Airborne Laser Scanning data (OPALS) program [50]. The ALS layers were selected based on their quality: if pixels on vegetation were identified as non-values or a "salt and pepper" effect was visible, the products were not further analysed. The "salt and pepper" effect is the occurrence of noise in the image in the form of isolated pixels with high local spatial heterogeneity between adjacent pixels [51]. We selected four layers from the BCAL Intensity Products: maximum, minimum, mean, and standard deviation. From among the BCAL Vegetation Products, 12 layers related to the height and structure of vegetation cover were used: minimum, maximum, range, mean, and coefficient of variation, alongside seven height percentiles products (the 5th, 10th, 25th, 50th, 75th, 90th, and 95th percentiles of all height points). The OPALS attributes for vegetation were selected, including the normalised mean, median, minimum, maximum, variance, RMS, and quantiles: 0.05, 0.10, 0.25, 0.75, 0.90, and 0.95. For all points, we also chose 4 amplitude parameters: mean, maximum, minimum, and median.
Finally, the following datasets for each campaign were created: hyperspectral bands (HS-430 layers), remote sensing vegetation indices (VI-49 layers), BCAL (4 layers for Intensity and 12 layers for Vegetation), and OPALS (17 layers). All layers were saved with a 1 m resolution corresponding to the HySpex images.

Statistical Analysis
The next step was an analysis of the differences between R. caesius and other types of vegetation based on the calculated remote sensing layers. We also sought to assess the differences between R. caesius and other types of vegetation. Rubus caesius polygons were divided into three classes for each campaign based on a species percentage coverage of 40-60%, 70-80%, or 90-100% (Table 2). Other classes (non-Rubus) were divided based on dominant species (e.g., Cirsium arvense, Solidago spp., and Tanacetum sp.) or for polygons without clear one dominant species based on vegetation classes distinguished according to Matuszkiewicz [17] (Artemisietea, Molinio-Arrhenatheretea, Phragmitetea). The above-mentioned round reference vector polygons were rasterised based on remote sensing layers, and pixels that fit entirely within the polygons were used in further analyses. The analyses were performed by comparing classes of R. caesius (coverage 40-60%, 70-80%, and 90-100%) to other vegetation classes (non-Rubus). In the early summer campaign, three R. caesius classes were analysed, whereas in the summer and autumn campaigns, only two were analysed: coverage 70-80% and 90-100% (there were not enough polygons for the R. caesius 40-60% class). Based on the number of pixels, six different vegetation classes were chosen for analysis: Artemisietea vulgaris, Molinio-Arrhenatheretea, Phragmitetea, Cirsium arvense, Solidago with Tanacetum, and trees.
To find the most differential layers, a two-stage process was performed for each of the three campaigns. First, all available layers in a given dataset (HS, VI, BCAL, or OPALS) were inspected for a correlation between their bands using a Tau-Kendal test. All layers with a pairwise correlation coefficient higher than 0.95 were removed. Next, Remote Sens. 2021, 13, 107 8 of 22 the remaining layers in each dataset were analysed using Linear Discriminant Analysis-LDA [38,40,52]. In the process, forward variable selection for classification using LDA was performed for each analysed dataset. Thus, for each dataset, we were able to rank layers based on their order of selection by LDA and their overall contribution to the correctness rate. The maximum number of layers selected by LDA was limited to 40. During the experiments, none of the datasets had more than 40 layers marked as important. Next, all layers deemed to be important were subject to a Non-Parametric Multivariate Analysis of Variance (NPMANOVA) to determine how good the separation was between classes in the given layers selected by LDA. All analyses were performed in the R software environment [53] using the packages caret [54], klaR [55], and vegan [56].
Different scenarios were analysed, including three R. caesius classes dependent on coverage compared to non-Rubus classes in three campaigns on four datasets: spectral reflectance (HS), vegetation indices (VI), OPALS, and BCAL. For each of the datasets, 100 iterations of the stratified random sampling procedure were performed, where 30 pixels from the R. caesius class and 30 from each of the non-Rubus classes were subsampled for decorrelation, LDA layer selection, and NPMANOVA. The pixels were randomised independently of the reference polygons so that each pixel could come from a different polygon or several pixels could be drawn within one polygon. The number of pixels was limited to 30 due to the number of pixels in each class. After analysis, each dataset was characterised by three values: the LDA correctness rate, the NPMANOVA F-value, and the layers that were selected by LDA to be the most useful. Based on the 100 iterations (for each the values of correctness rate, f values, differentiating layers), we identified the most differentiating layers. After 100 iterations of the LDA and NPMANOVA tests the F-values, correctness rates, and useful layer lists were compiled and visualised.
The next step was to analyse the results. Further analyses and interpretations were focused on the best datasets based on the correctness rates and F-values. The values for the BCAL and OPALS layers were noticeably lower compared to the hyperspectral bands and vegetation indices. For 430 spectral bands, the ranges of the spectrum were analysed for each campaign. This study assumed that the wavelength was useful for differentiating three R. caesius classes from non-Rubus if the band was within 10% of the highest frequency of occurrence. These ranges of the spectrum were analysed.
For the VI dataset, we analysed how many of the 100 iterations of each vegetation index were selected as differentiating. The same applied to HS, for which 10% of the most differentiating bands were analysed. Next, for the selected indices, maps of spatial distribution and values in each campaign were created. Based on this information and botanical knowledge, the data were analysed and interpreted.

Hyperspectral and ALS Data Comparison
The tests indicated the datasets and layers within them that differentiate R. caesius classes from non-Rubus classes. Based on the correctness rate from LDA and the F-value from NPMANOVA, the values were noticeably lower for the ALS data (the BCAL and OPALS datasets) ( Table 3 and Figure 4).
The correctness rate each time for AS products was equal to around 0.857, and there were no differences between iterations and different coverage levels of R. caesius. Even though the correctness rate values for ALS were high, this rate was always lower than that for the HS and VI. Based on this information, two datasets, BCAL and OPALS, were excluded from further analysis.
Higher values of the correction rate and F-values and a higher spread were acquired for the hyperspectral data (HS and VI). This differentiation was the greatest for the earlysummer campaign (C1), regardless of the coverage class, whereas the smallest was found for the autumn campaign (C3). However, the differences did not exceed a few percent. Higher values were calculated for classes with higher coverage, and the values were also less diverse. The average number of layers used to effectively differentiate the R. caesius Remote Sens. 2021, 13, 107 9 of 22 class from other vegetation for 100 iterations was similar for all datasets-from 10 to 14. Moreover, for the VI dataset, the F-value was generally higher compared to that of HS.

Spectral Band Dataset
As a result of the LDA analysis, we determined which spectral bands were the most useful for differentiating R. caesius depending on the species coverage and data acquisition date. Considering the spectrum bands with the highest frequency of occurrence as differentiating (10% of all bands), it can be concluded that the spectrum ranges are similar for the three campaigns and varied coverage (Figures 5 and 6).

Vegetation Index Dataset
For the VI dataset, we selected the indices that differentiate R. caesius classes the most based on the frequency of occurrence in each campaign (Figure 7). There were no significant differences between coverage classes and campaigns. The analysis showed that Normalised Difference Nitrogen Index (NDNI) has the highest potential to differentiate R. caesius from other types of vegetation among all analysed indices; in each campaign, at least once, the value of the frequency was above 90. High values were observed for Atmospherically Resistant Vegetation Index (ARVI) and Green Difference Vegetation Index (GDVI), which describe vegetation conditions, and for Cellulose Absorption Index (CAI), which identifies cellulose content. The indices that identify pigment content also had high frequency in each campaign, i.e., Anthocyanin Reflectance Index 1 (ARI1) and Caroten Reflectance Index 1 (CRI1). Figure 6. Average spectrum for different classes acquired from HySpex images and the bands (marked in dark grey, medium grey, and light grey dependent on Rubus coverage) that were defined as differentiating based on LDA analysis. The spectrum range was defined as differentiating when the frequency of the occurrence was above 15 (10% of 100 iterations). Average spectrum for different classes acquired from HySpex images and the bands (marked in dark grey, medium grey, and light grey dependent on Rubus coverage) that were defined as differentiating based on LDA analysis. The spectrum range was defined as differentiating when the frequency of the occurrence was above 15 (10% of 100 iterations).
The differentiating ranges are in VNIR (400-1000 nm) and SWIR (1000-2500 nm). Based on the values of frequency of occurrence, we located bands that were most commonly defined as differentiating in VNIR: There were only slight differences between the campaigns. For the early summer campaign (C1, the ranges of 920-970 nm and 1339 nm were used; for the summer campaign (C2), ranges of 774-777 nm were used; and for the autumn campaign (C3), we used 1720 nm, which is the absorption band for cellulose and lignin. For the early summer (C1) and summer (C2) campaigns, SWIR bands around 1665 nm and in the range 1653-1659 nm were observed to differentiate.

Vegetation Index Dataset
For the VI dataset, we selected the indices that differentiate R. caesius classes the most based on the frequency of occurrence in each campaign (Figure 7). There were no significant differences between coverage classes and campaigns. The analysis showed that Normalised Difference Nitrogen Index (NDNI) has the highest potential to differentiate R. caesius from other types of vegetation among all analysed indices; in each campaign, at least once, the value of the frequency was above 90. High values were observed for Atmospherically Resistant Vegetation Index (ARVI) and Green Difference Vegetation Index (GDVI), which describe vegetation conditions, and for Cellulose Absorption Index (CAI), which identifies cellulose content. The indices that identify pigment content also had high frequency in each campaign, i.e., Anthocyanin Reflectance Index 1 (ARI1) and Carotenoid Reflectance Index 1 (CRI1). When analysing the 10% most frequent vegetation indices in each campaign and R. caesius coverage with the highest frequency of occurrence, seven indices were noted: ARI1, CRI1, ARVI, GDVI, CAI, NDNI, and Modified Red Edge Simple Ratio Index (MRESR). However, the last index was observed only for R. caesius coverage of 70-80% for C3 (Table  4). Based on this information, six indices apart from MRESR were chosen for further comparisons (Figure 8). When analysing the 10% most frequent vegetation indices in each campaign and R. caesius coverage with the highest frequency of occurrence, seven indices were noted: ARI1, CRI1, ARVI, GDVI, CAI, NDNI, and Modified Red Edge Simple Ratio Index (MRESR). However, the last index was observed only for R. caesius coverage of 70-80% for C3 (Table 4). Based on this information, six indices apart from MRESR were chosen for further comparisons (Figure 8). Table 4. Vegetation indices (with information about spectrum used to calculate them) listed and marked with "+" among the 10% most differentiating layers for each campaign (C1, C2, C3) and specified Rubus caesius coverage.  The information on vegetation conditions was based on two analysed indices, ARVI and GDVI, which showed similar tendencies (Figures 8 and 9). The ARVI for R. caesius and most other tested types of vegetation had the highest values in early summer (C1); the average for all classes was 0.87, and the lowest in autumn (C3) was 0.73. Regardless of the examined term, the ARVI reached higher values for R. caesius (regardless of the percentage coverage) compared to other types of vegetation (except for patches with Solidago with Tanacetum) in C2. Within individual campaigns, the plots with varied R. caesius percentages differed only by 0.01. The GDVIs for R. caesius and for other types of vegetation were maximum in early summer (C1) and minimum in autumn (C3) (Figure 8 and Figure S2). The range of variation of this index for R. caesius during the growing season was 3308 points, from 6187 in C1 to 2879 in C3. Importantly, this index reached, on average, higher values for R. caesius (regardless of coverage) than for other types of vegetation. Only in early summer was the GDVI value for the vegetation of Artemisietea vulgaris class similar to that of R. caesius, while in autumn, only the value for rushes (Phragmitetea) was similar.

VI
The GDVI values for the three R. caesius coverage classes showed the same tendencies and slight differentiation within the campaign (maximum 190 points for C1). The indices related to pigment content generally had the highest values in summer. The ARI1 for the R. caesius class described the highest values in summer (C2), with 0.000579 at 70-80% coverage, and the lowest values in early summer, with 0.000314 for 40-60% coverage (Figures 8 and S1). In each of the three campaigns, we observed the highest values of the index for coverage of 70-80%. In early summer, for the R. caesius classes, this index had slightly higher values than those of non-Rubus types of vegetation. In summer and autumn (C2 and C3), the values were different in relation to the other types of non-forest vegetation, but differences were also observed for the three R. caesius classes. Only trees were characterised by a significantly higher value of ARI1. The CRI1 had the highest values for the R. caesius class in summer (C2), with 0.0022, and the lowest values in autumn (C3), with 0.0015 (Figure 8 and S4). For all tested types of vegetation, The indices related to pigment content generally had the highest values in summer. The ARI1 for the R. caesius class described the highest values in summer (C2), with 0.000579 at 70-80% coverage, and the lowest values in early summer, with 0.000314 for 40-60% coverage (Figure 8 and Figure S1). In each of the three campaigns, we observed the highest values of the index for coverage of 70-80%. In early summer, for the R. caesius classes, this index had slightly higher values than those of non-Rubus types of vegetation. In summer and autumn (C2 and C3), the values were different in relation to the other types of non-forest vegetation, but differences were also observed for the three R. caesius classes. Only trees were characterised by a significantly higher value of ARI1. The CRI1 had the highest values for the R. caesius class in summer (C2), with 0.0022, and the lowest values in autumn (C3), with 0.0015 (Figure 8 and Figure S4). For all tested types of vegetation, this index provided maximum values in summer. There was also a significantly greater variability of the index for the coverage class of 40-60% R. caesius than for higher coverage classes. In autumn, this species with a coverage of 70 to 100% was distinguished by a significantly higher value of CRI1 than most other types of examined non-forest vegetation. Only patches with Solidago and/or Tanacetum reached higher values in C3.
The values of NDNI, which shows the nitrogen content, were characterised by low variability within and between campaigns and presented values below 0 (Figure 8 and Figure S5). The values of this index ranged from −0.0381 for R. caesius in campaign C1 to -0.0647 for trees during campaign C2 (summer). The lowest value of this index was found in summer (C2). The highest values, depending on the type of vegetation, were found in early summer or autumn. The biggest differences between Rubus caesius and other vegetation classes were found in the summer campaign. For R. caesius, with coverage of 70-100%, the values were the highest apart from Molinio-Arrhenatheretea.
The cellulose content based on CAI, regardless of the campaign, was, on average, lower for R. caesius than for other types of vegetation ( Figure 8 and Figure S4). The values of this index were the lowest in autumn (C3), with an average 14 points, and the highest in summer (C2), with an average of 25 points. The difference between the values for R. caesius and other types of vegetation was the highest in early summer (C1). The value of the CAI index decreased with R. caesius coverage in each pixel.

Remote Sensing Data that Differentiate Rubus Caesius from Its Background
The performed analysis demonstrated that HS data and their products are useful for differentiating R. caesius from other vegetation types, whereas the results for ALS data were good but less useful in this application (Table 3 and Figure 4). Based on spectral information, it is possible to describe the biophysical components of the plants; the composition differs between species, so it is possible to use this kind of data for plant phenotyping [57]. Rubus is an expansive species that has different values for its biophysical variables (e.g., pigment content) and different vegetation condition (phenology) compared to other vegetation. This is why the VI, which indicates the status of the plants, differentiated R. caesius from non-Rubus vegetation.
The results of the LDA and NPMANOVA tests indicate that it is possible to differentiate the classes based on ALS products into BCAL and OPALS, but it is much more effective to use the HS bands or VI. Similar conclusions were obtained for Spiraea tomentosa mapping using ALS and hyperspectral data [58]. The ALS data were used to model the 3D construction of the vegetation canopy. Rubus caesius is surrounded by vegetation with a similar structure, so differentiation between species is possible only to a certain extent on the basis of parameters such as vegetation height and related derivatives. In this study, ALS data were used with a density of 7 points/m 2 , and all products were calculated with a 1-m spatial resolution. Using a higher point density or changing other ALS parameters could be useful for acquiring more detailed information. Moreover, different ALS products could be more successful than those used in this study. The ALS products were used for Rubus identification combined with HS data [32][33][34]. ALS data have rarely been used for the classification of non-forest vegetation. Based on ALS derivatives, were classified 24 types of vegetation in the Natura 2000 habitats, including salt marshes and steppic grasslands [59,60]. In the analyses, we used 18 variables from an ALS acquired in the growing season and three variables based on a DTM acquired in the leaf-off season. One of the key factors of effective differentiation was height. At the same time, hyperspectral data and ALS products were tested for detecting Spiraea tomentosa [58]. Using only the ALS product classification resulted in low accuracy (F1 for the S. tomentosa equalled about 42% depending on the time of data acquisition), while using only Minimum Noise Fraction bands, the value of F1 was at least 77%. A combination of these data resulted in a lower accuracy than that based only on the MNF bands (F1 about 68%).
In summary, for Rubus differentiation, ALS parameters are less useful compared to optical data probably because the height of Rubus is similar to that of the surrounding vegetation. However, based on the previous research, it can be stated that the combination of the two types of data increases the accuracy of the classification. Data fusion was tested for R. caesius classification [31] as well as for different types of vegetation [61][62][63]. At the same time, the combination of hyperspectral and ALS data was not tested in this study. The use of different types of data generally increases the cost and time of classification, which makes the monitoring more difficult.
The differentiation results in this study were slightly better for early summer and became worse with the flowing growing season. However, the differences were not significant, so it can be assumed that it is possible to distinguish Rubus from other species during most of the growing season. The previous results indicated summer to be the best time for data acquisition [35]. The classification of Rubus cuneifolius was performed in four seasons based on Landsat 8 and Sentinel-2 data. The obtained accuracies were not high, but there were some differences between the campaigns; the best results were acquired in the summer, and the worst were obtained in the winter.

Relation between Remote Sensing Data and the Functional Traits of Rubus Caesius
Based on the performed analysis, the results were found to be repetitive in each campaign-the generally determined VIs were based on the spectrum ranges defined as differentiating. There were slight shifts in the spectrum between campaigns, but most of the significant bands were similar. The shifts may be due to changes in incoming radiation or, less likely, to imperfections in the atmospheric correction. Two out of the seven most differentiating indices, ARVI and GDVI (Figure 8), used wide spectral ranges in the visible (red and blue) and near infrared spectra and were also indicated as discriminating based on LDA ( Figure 6). The wavelengths that were used to calculate the VI of the pigment content (ARI1-700 nm and CRI1-510 nm) and used in the CAI (2000, 2200 nm) and NDNI (1510, 680 nm) indices were partially defined as differentiating.
Differences were observed between the conclusions from the HS and VI datasets for ranges of 1968-2040 nm and bands above 2300 nm. The water content indices were not highlighted as the most differentiated indices. Water absorption bands are located in this spectrum range, so this range may be important for differentiating R. caesius from non-Rubus; however, these bands also had more noise than others.
A comparison between the three campaigns highlighted the potential weaknesses of the method. For the early summer campaign (C1), the significant ranges were around 960 nm. At the same time, these bands had significant noise due to errors in data acquisition, which is visible in the spectral profile ( Figure 6). The errors are probably related to the gap in spectrum during image acquisition and therefore, difficulty in combining data from two scanners (HySpex VNIR-1800 and SWIR-384). It can be concluded that in the case of acquisition errors, such a range could be incorrectly indicated as differentiating.
Based on the acquired results, it can be stated that four out of the seven most differentiating VIs had indicative values for R. caesius during the whole vegetation period: ARVI, CAI, CRI1, and NDNI ( Table 4). The other three were only indicative in individual campaigns, with ARI1 in early summer (C1), GDVI in summer (C2), and MRESR in autumn (C3).
The highest information potential was presented by the NDNI, which assumes average higher values for R. caesius than for the non-Rubus vegetation throughout the growing season ( Figure 8). The NDNI shows the nitrogen content in the leaves, which was proven for species growing in natural ecosystems such as the common reed (Phragmites australis) [64] and crops [65]. The large separability potential for the NDNI is likely due to the fact that R. caesius is a nitrophilous species [66] and may have a higher nitrogen content in its leaves than non-Rubus species. The maximum average NDNI value was additionally observed in early summer when, according to the results of the previous research, Rubus has the highest nitrogen content in tissues [67].
Another VI with a high separability potential is ARVI, which describes the plant's condition when reducing the influence of the atmosphere. This VI was used, among others, to identify Xylella fastidiosa infection in olive orchards [68]. On average, reference polygons with R. caesius have higher values of ARVI then non-Rubus species throughout the growing season and decrease from C1 to C3. High values of this VI are associated with high chlorophyll content in dewberry leaves [69]. CAI also has a high potential to identify Rubus, with average values significantly lower for Rubus than for the background (Figure 8). This lower value is related to the biology of the investigated species, where the CAI value is strongly related to the presence/visibility of the non-photosynthetic parts of plants [70], litter/dry biomass [69], and senesced vegetation [71]. The non-Rubus class consists mainly of non-forest communities constructed with perennials or annual species characterised by a relatively large coverage of dry biomass. On the other hand, the patches with R. caesius are characterised by high vigour and large coverage with horizontally oriented photosynthetic leaves. This was also confirmed by the high GDVI values compared to the background, which indicated a correlation with the Leaf Area Index [72]. Another important factor for distinguishing R. caesius from the background is the pigment indices such as CRI1, which describes the carotenoid content in plant tissues. The values of the CRI1 indicate that R. caesius is also characterised by a higher content of carotenoids compared to non-Rubus vegetation, which is related to the presence of high chlorophyll content, mainly to protect the photosynthetic system against photooxidation [69]. High carotenoid content is particularly visible in summer (C2) and autumn (C3) (Figure 8) for leaves exposed to full sun.

Remote Sensing Data Useful for Rubus Caesius Identification
Considering the aforementioned information, it can be assumed that the use of data with high spectral resolution is unnecessary because it is possible to detect R. caesius using only a few spectral ranges-much wider than those used in HS data. Based on the acquired results, it is possible to determine which ranges and bands allow one to differentiate R. caesius from non-Rubus. Considering the possibility of calculating the VI (ARVI and GDVI), it is important to use data in the visible and near infrared range. A red-edge band and potentially the water absorption range would also be useful.
Most of the previous studies on genus Rubus mapping were based on HS data. The following species were analysed with these data: R. armeniacus, R. fruticosus, and R. caesius [30][31][32][33][34]. This identification was mostly successful but also time consuming. It was possible to classify R. fruticosus L. agg. using HyMap [30] and Rubus spp. using HySpex [31]. Further, Rubus armeniacus was identified based on CASI and ALS products [32][33][34]. Again, these approaches were successful but also time-consuming, so their use in monitoring would likely be ineffective. Good results were acquired for Rubus cuneifolius (the best results among all tested methods were an OA of 82% and a kappa of 0.71) based on SPOT 5 [36]. In this research, the authors also suggested that successful mapping of the species is dependent, e.g., on the spatial and spectral resolution, patch size, time of acquisition, and other vegetation in the background.
An example of a sensor based on the aforementioned conditions that is not suitable for R. cuneifolius identification and has already been tested is Landsat (PA 52% and UA 45% for the best classification) [35]. The authors concluded that the images' spatial resolution was insufficient but also that the red-edge is not missing.
An example of a sensor that meets the aforementioned conditions is Sentinel-2, which offers a sufficient spectral resolution. The sensor also has red-edge and short-wave infrared bands. The spatial resolution makes it possible to identify R. caesius patches. The temporal resolution is also sufficient to acquire cloudless images during the whole growing season, including summer. Moreover, these data can be used for easy vegetation monitoring, as they are provided free of charge. Using this sensor, R. cuneifolius was classified with a PA of 80% and a UA of 45% [35]. The authors proved that the Sentinel-2 bands of NIR, red-edge, and SWIR are crucial in mapping, which is also a conclusion in the present study. Higher accuracy could possibly be achieved by using vegetation indices calculated based on Sentinel-2. However, the resolution of the Sentinel-2 data may be a limitation to its use. In case of dominance of small patches (below 20 m 2 ) linear spectral unmixing might be necessary to identify the Rubus. WorldView-3 could be a useful sensor for identification with a higher spatial resolution and at the same time potentially sufficient spectral resolution. It was already successfully used to with GEOBIA to extract olive tree crowns [73], oil palms [74], or broadleaf trees [75].

Conclusions
Based on the performed analysis, Rubus caesius can be distinguished from the background using optical data. Differentiation based on LDA and NPMANOVA analyses were successful for each of the analysed campaigns in early summer, summer, and autumn. ALS data were less useful for identification, which may be related to the insufficient variation between Rubus and non-Rubus classes in vegetation height.
Differentiation was possible for each R. caesius coverage class (40-60%, 70-80%, and 90-100%), but higher values of the correctness rate were observed for coverage of 90-100%. It was proven that regardless of the time of image acquisition, the selected spectral ranges were sufficient to differentiate R. caesius from non-Rubus. The analysis indicated what spectral ranges can be useful: visible (especially blue and red light), red-edge, near infrared, and possibly SWIR ranges are essential. It was proven that certain vegetation indices can differentiate R. caesius from other species, including ARI1, CRI1, ARVI, GDVI, CAI, NDNI, and MRESR. Only NDNI differentiated species regardless of the acquisition date and R. caesius coverage, while the other indices were not useful for every analysed scenario. For example, ARI1 could be used in early summer, and GDVI could be used in summer. The VI values were also found to be consistent with Rubus' biophysical properties. On this basis, it can be concluded that the selected variables were not random and can be used in classifications for different areas and times.
The obtained results indicate that it is possible to classify R. caesius using images with a lower spectral resolution than HS data. The high spectral and radiometric resolution makes it possible to identify brambles with high accuracy, but a great deal of information must be processed using dedicated algorithms, which takes considerable time. At the same time, identification of the species in monitoring must be precise and as fast and easy as possible. Moreover, the procedure should be repeatable in different areas and as simple as possible to carry out. Monitoring protected vegetation areas requires easy species identification, especially identification of species that pose a threat to natural habitats and valuable plant species. The spatial extent of R. caesius as an expansive species, should be monitored, as it poses a serious threat to biodiversity.
There is a need for interdisciplinary collaborations between ecologists and RS specialists to develop new methods for studying invasion [76]. Combining remote sensing technology with botanical knowledge will help us understand the interactions between radiation and vegetation and to map functional traits. This will enable information to be scaled over large areas.
The next steps for researching Rubus caesius mapping should be performed based on multispectral images, which is consistent with the specified assumptions. It should also be tested whether these data will be useful in other areas, with different sizes of patches, or for other Rubus species.
Supplementary Materials: The following are available online at https://www.mdpi.com/2072-4 292/13/1/107/s1, Table S1: Remote sensing vegetation indices used in the study. R-reflectance value; Figure S1: The spatial distribution of ARI2 classes for three campaigns in the selected area; Figure S2: The spatial distribution of ARVI classes for three campaigns in the selected area; Figure  S3: The spatial distribution of CAI classes for three campaigns in the selected area; Figure S4: The spatial distribution of CRI1 classes for three campaigns in the selected area; Figure S5: The spatial distribution of NDNI classes for three campaigns in the selected area.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to restrictions of founding institution NCBR.