Phenology-Based Mapping of an Alien Invasive Species Using Time Series of Multispectral Satellite Data: A Case-Study with Glossy Buckthorn in Qu é bec, Canada

: Glossy buckthorn ( Frangula alnus Mill.) is an alien species in Canada that is invading many forested areas. Glossy buckthorn has impacts on the biodiversity and productivity of invaded forests. Currently, we do not know much about the species’ ecology and no thorough study of its distribution in temperate forests has been performed yet. As is often the case with invasive plant species, the phenology of glossy buckthorn di ﬀ ers from that of other indigenous plant species found in invaded communities. In the forests of eastern Canada, the main phenological di ﬀ erence is a delay in the shedding of glossy buckthorn leaves, which occurs later in the fall than for other indigenous tree species found in that region. Therefore, our objective was to use that phenological characteristic to map the spatial distribution of glossy buckthorn over a portion of southern Qu é bec, Canada, using remote sensing-based approaches. We achieved this by applying a linear temporal unmixing model to a time series of the normalized di ﬀ erence vegetation index (NDVI) derived from Landsat 8 Operational Land Imager (OLI) images to create a map of the probability of the occurrence of glossy buckthorn for the study area. The map resulting from the temporal unmixing model shows an agreement of 69% with ﬁeld estimates of glossy buckthorn occurrence measured in 121 plots distributed over the study area. Glossy buckthorn mapping accuracy was limited by evergreen species and by the spectral and spatial resolution of the Landsat 8 OLI. in glossy buckthorn cover, forest cover types, site conditions (e.g., soil drainage, surficial deposit), and stand attributes. Site conditions and stand attributes (e.g., cover type, crown cover) were extracted from a 1:20 k provincial forest map [37]. We also aimed at maximizing the spatial distribution of the plots over the study area. We anticipated that evergreen foliage would reduce our capacity to isolate glossy buckthorn phenological characteristics from those of coniferous trees. For this reason, we tried to limit the number of plots in forest stands dominated by coniferous tree species. Finally, all located at least each other or from any


Introduction
Biological invasions by alien species can alter biodiversity and ecosystem functioning [1], with potential negative consequences on ecosystem services as well as on agricultural and forest resources [2]. The impacts of alien plants are observed at the species, community, and ecosystem levels [3]; for

Field Sampling
During the summer of 2015, we assessed glossy buckthorn cover in 124 sample plots (30 × 30 m) distributed over 45 private forest lots of the Richmond and Cookshire zones. A selection of candidate forest stands for field sampling was made based on information on glossy buckthorn history and cover obtained from the agency responsible for managing private forest land in the study region. Two criteria were used for identifying candidate stands: firstly, glossy buckthorn had to cover at least 900 m 2 of the stand total area and, secondly, observations of its presence in the stand needed to go back to at least three years prior to field sampling. These criteria ensured coherence between field measurements and Landsat 8 images and also increased the likelihood of detecting glossy buckthorn in 30 × 30 m pixels. From the set of candidate stands, we made the final selection with the objective of sampling the widest possible range in glossy buckthorn cover, forest cover types, site conditions (e.g., soil drainage, surficial deposit), and stand attributes. Site conditions and stand attributes (e.g., cover type, crown cover) were extracted from a 1:20 k provincial forest map [37]. We also aimed at maximizing the spatial distribution of the plots over the study area. We anticipated that evergreen foliage would reduce our capacity to isolate glossy buckthorn phenological characteristics from those of coniferous trees. For this reason, we tried to limit the number of plots in forest stands dominated by coniferous tree species. Finally, all plots were located at least 50 m from each other or from any

Field Sampling
During the summer of 2015, we assessed glossy buckthorn cover in 124 sample plots (30 × 30 m) distributed over 45 private forest lots of the Richmond and Cookshire zones. A selection of candidate forest stands for field sampling was made based on information on glossy buckthorn history and cover obtained from the agency responsible for managing private forest land in the study region. Two criteria were used for identifying candidate stands: firstly, glossy buckthorn had to cover at least 900 m 2 of the stand total area and, secondly, observations of its presence in the stand needed to go back to at least three years prior to field sampling. These criteria ensured coherence between field measurements and Landsat 8 images and also increased the likelihood of detecting glossy buckthorn in 30 × 30 m pixels. From the set of candidate stands, we made the final selection with the objective of sampling the widest possible range in glossy buckthorn cover, forest cover types, site conditions (e.g., soil drainage, surficial deposit), and stand attributes. Site conditions and stand attributes (e.g., cover type, crown cover) were extracted from a 1:20 k provincial forest map [37]. We also aimed at maximizing the spatial distribution of the plots over the study area. We anticipated that evergreen foliage would reduce our capacity to isolate glossy buckthorn phenological characteristics from those of coniferous trees. For this reason, we tried to limit the number of plots in forest stands dominated by coniferous tree species. Finally, all Remote Sens. 2020, 12, 922 4 of 15 plots were located at least 50 m from each other or from any spectrally contrasting feature, such as roads and clear-cuts, to reduce potential edge effects on OLI reflectance values [38].
We assessed plot-level glossy buckthorn cover using a modified version of the point-intercept method [39]. Plots were measured between May 12 and June 6 of 2015. In each plot, we established three 30-m long transects that were spaced 10 m apart and oriented north-south. Along each transect, we established one 4 m 2 -circular subplot (diameter = 1.13 m) every 5 m (7 subplots/transect), in which we assessed the occurrence of glossy buckthorn (presence/absence) and of other tree or shrub species in the subplot. Plot-level glossy buckthorn cover (%) was then calculated as the ratio of the number of subplots containing glossy buckthorn to the total number of subplots (N = 21). All sampled plots were geolocated using a Garmin Oregon 650 GPS (the Horizontal Dilution of Precision (HDOP) under the forest canopy was below 6 m using the Wide Area Augmentation System, or WAAS).

Processing of Landsat-8 OLI Images
We downloaded a time series of Landsat 8 OLI images (L1T) encompassing our study area from the USGS Global Visualization Viewer (GloVis) website (https://glovis.usgs.gov/). We limited our time series to images acquired between April and November during the years 2013 to 2015 and rejected images with a cloud cover above 10% (Table 1). Differences in atmospheric conditions between images in the time series can reduce the performance of change detection algorithms [40,41]. To minimize this source of noise, we used the ATCOR Ground Reflectance module in PCI Geomatica [42] to transform OLI digital numbers (DN) into surface reflectance values. Standard continental atmospheric properties and rural aerosol conditions were assumed, and vegetation was considered as a Lambertian surface. Topographic effects were also corrected in ATCOR using the Canadian Digital Elevation Model [43]. Finally, we extracted surface reflectance values in OLI 30 m bands 4 (red, 0.64-0.67 µm) and 5 (near-infrared, 0.85-0.88 µm) for all forest pixels in the time series. We masked non-forested pixels using the 1:20 k provincial forest maps.

NDVI Times Series
To identify the best temporal window for predicting glossy buckthorn occurrence and to obtain input and validation data for the phenology-based unmixing procedure, we calculated the NDVI for each date of OLI acquisition as: where ρ i is the surface reflectance in OLI band i. We generated NDVI time series for each measured plot by averaging, for each date, the NDVI values from a 3 × 3 pixel window centered on each plot coordinate. Finally, plot-level NDVI time series were used, along with glossy buckthorn percent cover measured in field plots, to determine the minimum cover value at which glossy buckthorn cover could be detected using the temporal unmixing approach.

Temporal Unmixing
The high level of spatial heterogeneity of forested areas, combined with the moderate spatial resolution of OLI, resulted in spectrally mixed forest pixels. To disentangle glossy buckthorn signals in mixed pixels, we used a linear spectral unmixing approach [44]. We applied spectral unmixing in the temporal domain, replacing spectral bands with image dates to retrieve the relative abundance (i.e., fraction) of cover types inside image pixels. Linear spectral mixture models consider pixel reflectance as a weighted sum of homogeneous spectral components called endmembers. Thus, we chose three temporally distinct endmembers (cover types), which also exhibit marked differences in their phenological trajectories: deciduous, coniferous, and glossy buckthorn. To create endmembers' temporal signatures, we used percent cover data measured in field plots. For each endmember, we selected one plot in which endmember cover type was at least 70% of the plot area. Field plots used for endmember selection were excluded from subsequent analyses. Endmember percentage cover values were 89%, 70%, and 80% for glossy buckthorn, deciduous, and coniferous endmembers, respectively. To minimize the influence of soil and background vegetation on endmember NDVI values, deciduous and coniferous endmember plots were selected among dense forest stands in which glossy buckthorn was absent. The glossy buckthorn endmember plot was located in a stand where tree cover was absent and that had the highest glossy buckthorn cover of all measured plots. To ensure a well-closed canopy, we used a summer NDVI image (August 19, 2015) to delineate, for each endmember, a group of contiguous pixels from the plot center. Endmember pixels were then extracted from each NDVI image and their descriptive statistics (minimum, maximum, mean, standard deviation) were extracted, by endmember and by date, to be used as inputs to the linear unmixing procedure. We used a t-test to confirm phenological differences between each endmember, which suggested that plots invaded by glossy buckthorn and those from which it was absent could be temporally distinguished. To account for the multiplicity of the test at each date (N = 3), we used the Bonferroni correction and used a threshold of 0.017 (0.05/3) to consider differences as significant (see Appendix A).
A linear spectral mixture model was created in the ENVI 4.8 software [45] and applied to unmix the time series of NDVI images. This model was used to assess the fraction of each endmember at the forested pixel level. The mixture model used in this procedure assumed that the NDVI time series results from each endmember. We used a fraction value of 0.4 to determine if a pixel was invaded by glossy buckthorn or not. This fraction was determined using a combination of the receiver operating characteristic (ROC) and Pareto front. The ROC curve was used to evaluate the effect of different fraction values on the false negative and true positive rates. We then used a Pareto front to determine the threshold fraction that allowed to minimize false negative rate while also maximizing the true positive rate (see Appendix B). The classification accuracy was evaluated using data from the remaining 121 field plots (109 invaded and 12 non-invaded). We used confusion matrices to calculate omission (lack of predictions) and commission (excess of predictions) errors and the level of agreement (accuracy) between field observations and model predictions.

Field-Based Measurements of Glossy Buckthorn Cover
The average glossy buckthorn cover from all field plots was 62.6%, with minimum and maximum values of 0% and 89%, respectively ( Table 2). Glossy buckthorn cover differed significantly between the two zones, with higher values in Cookshire than in Richmond (65.3% vs. 57.9%, respectively; p = 0.008).

Validation of the Phenology-Based Map of Glossy Buckthorn Occurrence
Glossy buckthorn was detected in only two of the 21 plots dominated by conifer trees, suggesting that dense conifer canopies may be an obstacle to the accurate detection of glossy buckthorn, regardless of the time of year. Therefore, coniferous-dominated plots were excluded, which left 100 plots from subsequent validation analyses. The confusion matrix (Table 3) shows that the classification method achieved a global accuracy of 69%, meaning that 69% of present glossy buckthorn was classified accurately (true positive rate = 62/90). Similarly, 70% of the sites in which glossy buckthorn was absent were correctly classified. On the other hand, our approach could not detect glossy buckthorn presence in 31% of the plots where it was present (false negative rate = 28/90), whereas in 30% of the cases its presence was predicted when it was actually absent from the plot (false positive rate = 3/10). Table 3. Confusion matrix for glossy buckthorn occurrence estimated from temporal linear unmixing of the Landsat-8 OLI normalized difference vegetation index (NDVI) image series. The threshold value (% cover) was used to determine if glossy buckthorn was present if the endmember fraction map was set to 0.4 (minimum glossy buckthorn cover measured value in presence plots was 27%).

Predicted
The important variability in glossy buckthorn percent cover among the 90 plots where it was found allowed splitting the dataset into three cover classes. The classification performance increased as the class lower bound increased, with an 11% increase from the class covering the largest range in buckthorn cover (≥25%) to that containing only plots with a high buckthorn cover (≥75%) ( Table 4). Table 4. Agreement and disagreement values calculated from the confusion matrices of different glossy buckthorn cover classes. Each confusion matrix used the same amount of absence plots (N = 10) but decreasing numbers of presence plots (≥ 25%: N = 90 (lowest glossy buckthorn cover), ≥65%: N = 60 (average glossy buckthorn cover), ≥ 75%: N = 32).

≥25%
≥65% ≥75% Given the differences in forest cover between the two zones (Table 2), we compared the performance of our classifier by calculating agreement/disagreement values separately for each zone. We found a marked difference in performance between Richmond (59% agreement) and Cookshire (75% agreement). In most forested areas of the Richmond zone, glossy buckthorn was identified by temporal unmixing analysis. There were some places where the probability of occurrence of glossy buckthorn was more than 90% ( Figure 3). In most forested areas of the Richmond zone, glossy buckthorn was identified by temporal unmixing analysis. There were some places where the probability of occurrence of glossy buckthorn was more than 90% (Figure 3).
The results show a different pattern for the Cookshire zone ( Figure 4). Indeed, glossy buckthorn was present but scattered in the landscape. There was some forested area where glossy buckthorn was likely absent. The results show a different pattern for the Cookshire zone ( Figure 4). Indeed, glossy buckthorn was present but scattered in the landscape. There was some forested area where glossy buckthorn was likely absent.

Discussion
An earlier study using Landsat-7 ETM+ data suggested that glossy buckthorn should compose at least 70% of the canopy in pixels to be successfully identified [46]. Likewise, we found that the mapping of buckthorn occurrence was more accurate in plots where glossy buckthorn cover was higher than 70%. High spatial resolution hyperspectral or multispectral images would have been a better tool than Landsat imagery for investigating zones where glossy buckthorn covered less than 70% of the canopy. However, the costs associated with high spatial resolution images (<1 m/pixel) are high and they generally require advanced processing methods and a high computing power [47]. Becker et al. [46] mapped glossy buckthorn and common buckthorn (Rhamnus cathartica L.) in a territory using 50 × 50 m experimental plots. Unfortunately, their classification with Landsat 5 TM and Landsat 7 ETM + was found to be ineffective in patches smaller than 50 × 50 m. Most of our field plots showed a dense 2-3-m high glossy buckthorn cover. We commonly encountered dispersed small patches of dense glossy buckthorn bushes, which were undoubtedly difficult to identify using 30-m Landsat 8 OLI images. It is likely that this particularity, combined with a well-developed tree cover, has contributed to the high false negative rate we obtained and that this would limit the precision of this approach if applied for mapping glossy buckthorn occurrence in stands where similar attributes exist. Our mapping was cautious because it included only places heavily invaded by glossy buckthorn and can be used as a reference map for future analysis. The high commission error rate reveals that the model failed to predict the absence of glossy buckthorn in the study area. Indeed, the model underestimated the actual invasion of glossy buckthorn in most coniferous stands. This underestimation is particularly important since conifers occupy 17% of the study area. Therefore, it appears that multispectral optical remote sensing data such as Landsat 8 OLI is not the optimal tool for mapping glossy buckthorn occurrence in these stands since these plants remain

Discussion
An earlier study using Landsat-7 ETM+ data suggested that glossy buckthorn should compose at least 70% of the canopy in pixels to be successfully identified [46]. Likewise, we found that the mapping of buckthorn occurrence was more accurate in plots where glossy buckthorn cover was higher than 70%. High spatial resolution hyperspectral or multispectral images would have been a better tool than Landsat imagery for investigating zones where glossy buckthorn covered less than 70% of the canopy. However, the costs associated with high spatial resolution images (<1 m/pixel) are high and they generally require advanced processing methods and a high computing power [47]. Becker et al. [46] mapped glossy buckthorn and common buckthorn (Rhamnus cathartica L.) in a territory using 50 × 50 m experimental plots. Unfortunately, their classification with Landsat 5 TM and Landsat 7 ETM + was found to be ineffective in patches smaller than 50 × 50 m. Most of our field plots showed a dense 2-3-m high glossy buckthorn cover. We commonly encountered dispersed small patches of dense glossy buckthorn bushes, which were undoubtedly difficult to identify using 30-m Landsat 8 OLI images. It is likely that this particularity, combined with a well-developed tree cover, has contributed to the high false negative rate we obtained and that this would limit the precision of this approach if applied for mapping glossy buckthorn occurrence in stands where similar attributes exist. Our mapping was cautious because it included only places heavily invaded by glossy buckthorn and can be used as a reference map for future analysis. The high commission error rate reveals that the model failed to predict the absence of glossy buckthorn in the study area. Indeed, the model underestimated the actual invasion of glossy buckthorn in most coniferous stands. This underestimation is particularly important since conifers occupy 17% of the study area. Therefore, it appears that multispectral optical remote sensing data such as Landsat 8 OLI is not the optimal tool for mapping glossy buckthorn occurrence in these stands since these plants remain invisible to the instrument throughout the year under dense conifer canopies. This limitation might become more important in the coming decades, as climate change will potentially favor the extension of glossy buckthorn distribution in the Canadian boreal forest (from where is it actually relatively absent).
There were several reasons why our mapping was more effective in Cookshire than in the Richmond zone. First, field observations by managers and foresters showed that the invasion was extensive in Cookshire. This may also be explained by the difference in stand cover of the validation sites, which was lower in Cookshire than in Richmond. For sure, it was more difficult to identify glossy buckthorn through a dense forest canopy, even knowing that the leaves have fallen at some point, due to the interference caused by the presence of branches and tree trunks. By comparison, a study using a time series of NDVI to perform the characterization of glossy buckthorn and common buckthorn, managed to successfully identify these two species in Ohio and Michigan (USA) with a kappa of 0.73 [46]. This seems greatly facilitated by the ecological characteristics of the study area, since communities were mainly oak grasslands, low-cover environments without conifer species. Our results suggest that phenology-based mapping using unmixing and the time series of multispectral satellite data may be more efficient to detect species at a lower threshold than traditional approaches.
Along with complications due to the tree canopy, there were also other interferences, such as the presence of exposed soil and of other plant species, such as ferns (primarily Onoclea sensibilis L.) or speckled alder (Alnus rugosa (Du Roi) Spreng.). Speckled alder is a shrub that grows in the understory, very similar to glossy buckthorn, also considered a companion species of the latter. The presence of this understory vegetation modified the reflectance signal of the canopy. According to some authors, these variations can be over 18% in the red band, and over 10% in the infrared band [48]. Some studies have even used this variability in the reflectance values to detect the type of understory vegetation [49][50][51]. Using Landsat 8 (OLI) allowed the classification of an invaded forest according to forest types, rather than the understory species themselves. The spatial resolution of Landsat 8 sensor (OLI) thus represents a limit when trying to map smaller entities than the pixel.
Other images with a higher spatial or spectral resolution (such as WorldView, IKONOS, Quickbird) could have been used in this study, but the costs associated with these images are prohibitive and, moreover, they may not be available for the desired dates and regions of interest. This is also the case for airborne (AVIRIS, CASI, etc.) or in-orbit (EO1-Hyperion) hyperspectral imagery. Hyperspectral sensors could greatly improve glossy buckthorn mapping because of their ability to discriminate pixels by taking advantage of the richness of the spectral information contained in large numbers of high-spectral resolution bands [52]. In addition, it is sometimes necessary to make several corrections to obtain images without blur, which complicates the use of these images [53]. Although they have proved to be efficient in mapping alien species, e.g., [54], the acquisition of these images by airborne sensors can be complex with respect to the availability of flying machines or current regulations [27]. LiDAR was used for classifying understory species in coniferous stands [55,56]. This tool could also improve outcomes for mapping glossy buckthorn, since the tool is able to provide information from under the tree canopy. Because the structure of coniferous stands differs from that of deciduous forests [57], and also because it can be modified by the presence of an alien invasive species in the understory [58], it is likely that the presence of glossy buckthorn in coniferous woods can be detected with LiDAR. In addition, the amount of commission errors in the absence class indicates that more invaded pixels were not properly classified as a presence, reflected again by an underestimation of the cover.

Conclusions
Temporal unmixing performed on a time series of Landsat8 (OLI) NDVI images allowed for the detection of glossy buckthorn in deciduous and mixed forest. The mapping allowed to identify forest stands heavily overgrown by glossy buckthorn, but did not allow detection in coniferous stands. Thus, the study helped to partially verify the hypothesis that phenology allows to discriminate the temporal signature of glossy buckthorn from that of other forest cover types found in the study area. In order to improve the occurrence mapping of glossy buckthorn, other avenues should be considered such as increasing the number of presence and absence sites to better reflect the great variability in forest structures and compositions (cover types, canopy closure, settlement patterns, etc.). This would allow for the selection of more representative points of absence for the calculation of homogeneous spectral components. Other spectral indices less likely to saturate in well-developed biomass stands (e.g., EVI or ARVI) could improve the glossy buckthorn detection threshold.  considered such as increasing the number of presence and absence sites to better reflect the great variability in forest structures and compositions (cover types, canopy closure, settlement patterns, etc.). This would allow for the selection of more representative points of absence for the calculation of homogeneous spectral components. Other spectral indices less likely to saturate in well-developed biomass stands (e.g., EVI or ARVI) could improve the glossy buckthorn detection threshold.  .  Table A1. Comparisons of endmember mean values: t-tests F-stats and P-values for each pair of endmember and the date of image acquisition. P-values in bold indicate a significant difference in mean between two groups for a given date. To account for the multiplicity of the test at each date (N = 3), we applied the Bonferroni correction. P-values below 0.017 (0.05/3) were considered as significant.  Appendix B Figure A2. The ROC curve used to determine the optimal threshold value for predicting the occurrence of glossy buckthorn. A Pareto front was used to determine the fraction of glossy buckthorn Figure A2. The ROC curve used to determine the optimal threshold value for predicting the occurrence of glossy buckthorn. A Pareto front was used to determine the fraction of glossy buckthorn that allows one to minimize the false negative rate and maximize the true positive rate. The black line describes the behavior of the ROC algorithm for all possible threshold values, while the blue dashed line represents a random classifier with no power of discrimination.