Green Vegetation Cover Dynamics in a Heterogeneous Grassland: Spectral Unmixing of Landsat Time Series from 1999 to 2014

: The ability to quantify green vegetation across space and over time is useful for studying grassland health and function and improving our understanding of the impact of land use and climate change on grasslands. Directly measuring the fraction of green vegetation cover is labor-intensive and thus only practical on relatively smaller experimental sites. Remote sensing vegetation indices, as a commonly-used method for large-area vegetation mapping, were found to produce inconsistent accuracies when mapping green vegetation in semi-arid grasslands, largely due to mixed pixels including both photosynthetic and non-photosynthetic material. The spectral mixture approach has the potential to map the fraction of green vegetation cover in a heterogeneous landscape, thanks to its ability to decompose a spectral signal from a mixed pixel into a set of fractional abundances. In this study, a time series of fractional green vegetation cover (FGVC) from 1999 to 2014 is estimated using the spectral mixture approach for a semi-arid mixed grassland, which represents a typical threatened, species-rich habitat in Central Canada. The shape of pixel clouds in each of the Landsat images is used to identify three major image endmembers (green vegetation, bare soil / litter, and water / shadow) for automated image spectral unmixing. The FGVC derived through the spectral mixture approach correlates highly with ﬁeld observations (R 2 = 0.86). Change in the FGVC over the study period was also mapped, and green vegetation in badlands and uplands is found to experience a slight increase, while vegetation in riparian zone shows a decrease. Only a small portion of the study area is undergoing signiﬁcant changes, which is likely attributable to climate variability, bison reintroduction, and wildﬁre. The results of this study suggest that the automated spectral unmixing approach is promising, and the time series of medium-resolution images is capable of identifying changes in green vegetation cover in semi-arid grasslands. Further research should investigate driving forces for areas undergoing signiﬁcant changes.


Introduction
Grasslands occupy nearly half of the terrestrial globe and play a critical role in supplying food, fiber, and fuel, supporting the biodiversity of animals and plants, maintaining water and air quality, and supporting ecological processes that sustain ecosystems and landscapes [1,2]. However, many grasslands are highly fragmented [3], and the remnants are vulnerable to the effects of adjacent land use, management regimes, and climate variability [4]. Protecting and monitoring grassland health has been a major priority for many local, national, and international agencies.
In summary, the improvement of radiative transferring modeling for mixed scenes, including both green and senescence vegetation, is well underway.
The above-mentioned approaches use an image pixel as the smallest entity, and the spatial heterogeneity within the pixel is ignored. An increase heterogeneity within mixed pixels not only reduces VI-regression accuracies [23] but also affects the precise simulation of the radiative characteristics of land surfaces and the inversion of land surface parameters in the radiative transfer models [32]. This is particularly true for the heterogeneous landscape of semi-arid grassland ecosystems (Figure 1), where patches are often smaller than the dimensions of a median spatial resolution image (e.g., Landsat 30 m) [33]. The relatively new method of the spectral mixture approach is a promising alternative and has improved to better characterize heterogeneous landscapes through the decomposition of a spectral signal from a mixed pixel into a set of fractional abundances, including FGVC using spectral signatures from pure ground components called endmembers. The spectral unmixing analysis describes physical components within the canopy and is thus different from the regression approach and appreciated for its generality in nature. Spectral unmixing may be applied to a wide range of landscapes over time, depending on the type of spectral unmixing algorithm employed. For semi-arid ecosystems, spectral unmixing of remote sensing data yielded satisfactory results. For example, Elmore et al. [13] quantified vegetation change in semi-arid environments and found that the spectral unmixing approach correctly determined 87% of the change in percent live vegetation cover, superior to NDVI, which only determined 67% of the change.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 20 radiative transferring modeling for mixed scenes, including both green and senescence vegetation, is well underway. The above-mentioned approaches use an image pixel as the smallest entity, and the spatial heterogeneity within the pixel is ignored. An increase heterogeneity within mixed pixels not only reduces VI-regression accuracies [23] but also affects the precise simulation of the radiative characteristics of land surfaces and the inversion of land surface parameters in the radiative transfer models [32]. This is particularly true for the heterogeneous landscape of semi-arid grassland ecosystems (Figure 1), where patches are often smaller than the dimensions of a median spatial resolution image (e.g., Landsat 30 m) [33]. The relatively new method of the spectral mixture approach is a promising alternative and has improved to better characterize heterogeneous landscapes through the decomposition of a spectral signal from a mixed pixel into a set of fractional abundances, including FGVC using spectral signatures from pure ground components called endmembers. The spectral unmixing analysis describes physical components within the canopy and is thus different from the regression approach and appreciated for its generality in nature. Spectral unmixing may be applied to a wide range of landscapes over time, depending on the type of spectral unmixing algorithm employed. For semi-arid ecosystems, spectral unmixing of remote sensing data yielded satisfactory results. For example, Elmore et al. [13] quantified vegetation change in semi-arid environments and found that the spectral unmixing approach correctly determined 87% of the change in percent live vegetation cover, superior to NDVI, which only determined 67% of the change. Figure 1. A typical mixed grassland scene with exposed soil, photosynthetic grass, and nonphotosynthetic grass.
The success of spectral unmixing is largely dependent on endmember selection. Earlier studies extracted endmembers from a spectra pool that was established by laboratory or fieldwork [34]; however, the spectral data in the pool are rarely acquired under the same conditions as the remote sensing data and may misrepresent real endmembers in the images. Recent studies have paid more attention to selecting quality endmembers from the images under investigation. These endmembers (i.e., image endmembers) address the limitation of utilizing lab or field-based spectra pool as they are collected under the same environmental conditions as mixed pixels in the images. It is widely noted that image endmembers should be selected from extreme values (i.e., the purest pixels and composed only for one member) of the image spectral feature space [35]. A typical mixed grassland scene with exposed soil, photosynthetic grass, and non-photosynthetic grass.
The success of spectral unmixing is largely dependent on endmember selection. Earlier studies extracted endmembers from a spectra pool that was established by laboratory or fieldwork [34]; however, the spectral data in the pool are rarely acquired under the same conditions as the remote sensing data and may misrepresent real endmembers in the images. Recent studies have paid more attention to selecting quality endmembers from the images under investigation. These endmembers (i.e., image endmembers) address the limitation of utilizing lab or field-based spectra pool as they are collected under the same environmental conditions as mixed pixels in the images. It is widely noted Remote Sens. 2020, 12, 3826 4 of 19 that image endmembers should be selected from extreme values (i.e., the purest pixels and composed only for one member) of the image spectral feature space [35].
Automated and fully objective methods for image endmember selection are essential for an efficient image unmixing practice. Several research groups have started to identify the image endmembers by investigating the shape of pixel clouds in spectral mixing spaces [36][37][38]. The optimal shapes of pixel clouds for spectral unmixing are identified to be triangles and tetrahedrons in 2-and 3-dimensional spectral mixing spaces, respectively [37]. When assigning the vertices of triangles and tetrahedrons as endmembers, spectral unmixing has been demonstrated to be valid and meaningful for all the pixels in the image. With this finding, the process of automated image spectral unmixing is thus narrowed down to the characterization of an appropriate spectral mixing space so that the shape of the pixel cloud resembles a triangle or a tetrahedron, and the vertices represent target endmembers.
This study aims to investigate the shape of pixel clouds in a series of Landsat images for a semi-arid grassland and use the identified image endmembers for automated image spectral unmixing. The ultimate goal is to provide an overview of the dynamics of fractional green vegetation cover for the semi-arid mixed grassland, which represents a typical example of threatened, species-rich habitat in central Canada.

Study Area
The study area was located in the West Block of the Grasslands National Park (GNP) (49 • 19 26.3" N, 107 • 32 37.4" W) in southwestern Saskatchewan, Canada ( Figure 2). This area is defined as a semi-arid, mixed-grass prairie ecosystem that falls within the northern extent of the Great Plains [7]. The park was established in 1988 to conserve and protect a portion of the remaining natural grasslands that were once so extensive in southwest Saskatchewan. The climate in the study area is typical of a semi-arid environment, with long cold and dry winters, short and hot summers, and low precipitation [39]. The mean annual temperature is approximately 4.1 • C and ranges from −10.8 • C in January to 18.5 • C in July, while the mean annual precipitation is approximately 352.5 mm. The soils in the West Block of the GNP are brown Chernozemic clay to clay loam soils according to the Canadian System of Soil Classification. The dominant native grass species found in the study area were June grass (Koeleria gracilis), needle-and-thread grass (Stipa comata), blue grama (Bouteloua gracilis), western wheatgrass (Agropyron smithii), and northern wheatgrass (Agropyron dasystachym). A non-native invasive grass species, crested wheatgrass (Agropyron cristatum), was also a dominant type found in the study area. Three typical grassland landscape units, including badlands, uplands, and riparian grasslands (Figure 2), can be found in GNP; however, the majority of the grassland landscape in GNP consists of exposed soil and mixed non-photosynthetic grass and green vegetation patches ( Figure 1).

Field Campaigns, Landsat Images, and Climatic Data
Seven field campaigns (i.e., 2003, 2005, 2006, 2010, 2011, 2012, and 2013) were carried out to determine various grassland biophysical characteristics during the peak growing season (Middle June to Middle July). A total of 118 field sites were visited over the seven years, and these sites were distributed across uplands, badlands, and riparian zone ( Figure 2). Samples in each site were collected along two 100 m transects that ran north-south and west-east, forming a perpendicular bisector. The perpendicular bisector design was developed for use with moderate spatial resolution sensor data such as Landsat, where each site will cover about nine pixels of a 30-m resolution Landsat image.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 20 Figure 2. Study area is located in the West Block of the Grasslands National Park, Saskatchewan, Canada. (a-c) are photos were taken from badlands, upland grassland, and riparian grassland areas, respectively.

Field Campaigns, Landsat Images, and Climatic Data
Seven field campaigns (i.e., 2003, 2005, 2006, 2010, 2011, 2012, and 2013) were carried out to determine various grassland biophysical characteristics during the peak growing season (Middle June to Middle July). A total of 118 field sites were visited over the seven years, and these sites were distributed across uplands, badlands, and riparian zone ( Figure 2). Samples in each site were collected along two 100 m transects that ran north-south and west-east, forming a perpendicular bisector. The perpendicular bisector design was developed for use with moderate spatial resolution sensor data such as Landsat, where each site will cover about nine pixels of a 30-m resolution Landsat image.
A detailed list of sites surveyed during each field campaign can be found in Table 1. To investigate the feasibility of Landsat images in constructing the spectral feature space and in distinguishing green vegetation from bare soil/litter and water for this grassland, we also collected a few field variables, including canopy composition and canopy hyperspectral reflectance. Field measurements were recorded within a 50 cm by 50 cm sampling frame along each 100 m transect at 10 m intervals, resulting in a total of 20 plots at each site. Canopy composition measurements, including percent coverage of the green vegetation, litter, soil, and water percentage coverages, were collected from each sampling site through both field interpretation and field photo classification. In specific, the field photos taken from each sampling frame were used to estimate the percent coverage of different features through unsupervised classification schemes in Geomatica ® 10.3 (PCI Geomatics, Richmond Hill, Ontario, Canada). The consistency between field interpreted vegetation cover, and the field photo classification was analyzed to ensure the validity of field data. In addition, field-based hyperspectral reflectance was measured from each plot to simulate Landsat bands and construct spectral feature space. The reflectance data were taken approximately 1 m above each plot using an Analytical Spectral Device (ASD) FieldSpec ® 3 Max field portable spectroradiometer (Malvern Panalytical Inc., Westborough, Massachusetts, USA), which is equipped with a fiber optic that has A detailed list of sites surveyed during each field campaign can be found in Table 1. To investigate the feasibility of Landsat images in constructing the spectral feature space and in distinguishing green vegetation from bare soil/litter and water for this grassland, we also collected a few field variables, including canopy composition and canopy hyperspectral reflectance. Field measurements were recorded within a 50 cm by 50 cm sampling frame along each 100 m transect at 10 m intervals, resulting in a total of 20 plots at each site. Canopy composition measurements, including percent coverage of the green vegetation, litter, soil, and water percentage coverages, were collected from each sampling site through both field interpretation and field photo classification. In specific, the field photos taken from each sampling frame were used to estimate the percent coverage of different features through unsupervised classification schemes in Geomatica ® 10.3 (PCI Geomatics, Richmond Hill, ON, Canada). The consistency between field interpreted vegetation cover, and the field photo classification was analyzed to ensure the validity of field data. In addition, field-based hyperspectral reflectance was measured from each plot to simulate Landsat bands and construct spectral feature space. The reflectance data were taken approximately 1 m above each plot using an Analytical Spectral Device (ASD) FieldSpec ® 3 Max field portable spectroradiometer (Malvern Panalytical Inc., Westborough, MA, USA), which is equipped with a fiber optic that has three detectors, each operating in either the visible, near-infrared (NIR), or short-wave infrared region. The ASD spectroradiometer was frequently calibrated with a certified white reference, and spectral measurements were only taken on days with clear sky conditions between 10:00 and 14:00 h. All field measured data were then amalgamated into mean values for each 100 by 100 m site, and the field green vegetation cover data for those sites were used to validate the Landsat-unmixed FGVC. Geometrically and atmospherically-corrected Landsat Level-2 images were downloaded between the years 1999 and 2014 from the USGS Landsat archive. Images with an acquisition time as close to peak growing season (15 June to 30 August) and as cloud-free as possible for the study area were eventually used in this work to minimize the effects of vegetation phenology and ensure consistency with field data (Table 1). This resulted in a time series of at least one image every summer, other than the year of 2004, in which year all images suffered from cloud issues, and as such no images were used from this year.
Daily precipitation data for this study were downloaded from a nearby weather station (Val Marie Southeast Saskatchewan, Climate ID: 4338412) from Environment Canada (http://climate.weather.gc. ca/historical_data/search_historic_data_e.html). The daily precipitation data for the weather station were accumulated to 3 days, 1 week, 2 weeks, 3 weeks, and a month before each image was taken to show how vegetation in badlands, uplands, and riparian zone responded to climate conditions.

Spectral Unmixing
The flowchart ( Figure 3) provides a general overview of the spectral unmixing procedure described below. Specifically, the spectral unmixing uses an image as input, followed by spectral mixing space construction, endmember identification, and linear unmixing. The resultant unmixing maps are then validated by ground truth data, and the final maps will then be produced along with an associated accuracy.

Determining an Appropriate Spectral Mixing Space
To determine which image bands can be used to distinguish photosynthetic grass from other typical ground features including litter, bare soil, and water, field measured spectra ( Figure 4a) were plotted against wavelengths and further used to simulate Landsat bands. Generally, photosynthetic grass has a typical green vegetation spectra curve: high reflectance in green, low in red, and very high in NIR. Litter and bare soil have a very similar spectra curve: reflectance showing an increasing trend from short to longer wavelength regions. In comparison with the spectra of typical green vegetation, the reflectance of a typical semi-arid mixed grassland canopy comprised of green vegetation, bare soil, and a large portion of senescent leaves is lower in green, higher in Red, and much lower in NIR. Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 20

Determining an Appropriate Spectral Mixing Space
To determine which image bands can be used to distinguish photosynthetic grass from other typical ground features including litter, bare soil, and water, field measured spectra ( Figure 4a) were plotted against wavelengths and further used to simulate Landsat bands. Generally, photosynthetic grass has a typical green vegetation spectra curve: high reflectance in green, low in red, and very high in NIR. Litter and bare soil have a very similar spectra curve: reflectance showing an increasing trend from short to longer wavelength regions. In comparison with the spectra of typical green vegetation, the reflectance of a typical semi-arid mixed grassland canopy comprised of green vegetation, bare soil, and a large portion of senescent leaves is lower in green, higher in Red, and much lower in NIR.
The spectral mixing space ( Figure 4b) constructed by simulated Landsat red and NIR bands from ASD spectroradiometer measurements indicates that these two bands are able to distinguish green vegetation from bare soil/litter and water. As identified in [40], the reflectance spectra of bare soil and litter are similar enough in simulated Landsat bands to be considered as one ground feature. The three ground features (green vegetation, bare soil/litter, and water) formed the vertices of a point cloud triangle in the red-NIR spectral mixing space ( Figure 4b). The majority of spectra from the mixed field plots are found to be within the triangle. The spectra in a few plots that are not within the triangle likely resulted from the different biophysical or/and biochemical properties of vegetation and soil within these plots, in comparison with those of the selected typical green vegetation and soil.  . Geometric interpretation of the mixture problem in a two-dimensional space using fieldcollected hyperspectral data; (a) field observed spectra curves of typical ground features including vegetation, litter, soil, water, and a mixed canopy; (b) point cloud plotted from field spectral simulated Landsat Thematic Mapper (TM) red and near-infrared (NIR) bands for the typical ground features (vertexes forming a triangle shade) and for all field-collected plot spectra (black dots).

Image Endmember Optimization and Spectral Mixing Space Translation
The above section demonstrated that simulated Landsat Red and NIR bands could be used to distinguish photosynthetic grass from other typical ground features, including litter, bare soil, and water. We can now proceed to identify endmembers from each image under study. The NIR and red values from all pixels in each Landsat image were plotted in the spectral scatterplot to form a specific spectral mixing space [36].
A recently proposed endmember optimization method was employed to automatically identify the endmembers of vegetation and bare soil/litter in the Landsat red-NIR feature space for each year [37]. Specifically, an endmember spatial measure index (ESM2) is directly proportional to the determinant (Det) of all endmember vectors, . Geometric interpretation of the mixture problem in a two-dimensional space using field-collected hyperspectral data; (a) field observed spectra curves of typical ground features including vegetation, litter, soil, water, and a mixed canopy; (b) point cloud plotted from field spectral simulated Landsat Thematic Mapper (TM) red and near-infrared (NIR) bands for the typical ground features (vertexes forming a triangle shade) and for all field-collected plot spectra (black dots).
The spectral mixing space (Figure 4b) constructed by simulated Landsat red and NIR bands from ASD spectroradiometer measurements indicates that these two bands are able to distinguish green vegetation from bare soil/litter and water. As identified in [40], the reflectance spectra of bare soil and litter are similar enough in simulated Landsat bands to be considered as one ground feature. The three ground features (green vegetation, bare soil/litter, and water) formed the vertices of a point cloud triangle in the red-NIR spectral mixing space (Figure 4b). The majority of spectra from the mixed field Remote Sens. 2020, 12, 3826 8 of 19 plots are found to be within the triangle. The spectra in a few plots that are not within the triangle likely resulted from the different biophysical or/and biochemical properties of vegetation and soil within these plots, in comparison with those of the selected typical green vegetation and soil.

Image Endmember Optimization and Spectral Mixing Space Translation
The above section demonstrated that simulated Landsat Red and NIR bands could be used to distinguish photosynthetic grass from other typical ground features, including litter, bare soil, and water. We can now proceed to identify endmembers from each image under study. The NIR and red values from all pixels in each Landsat image were plotted in the spectral scatterplot to form a specific spectral mixing space [36].
A recently proposed endmember optimization method was employed to automatically identify the endmembers of vegetation and bare soil/litter in the Landsat red-NIR feature space for each year [37]. Specifically, an endmember spatial measure index (ESM 2 ) is directly proportional to the determinant (Det) of all endmember vectors, where − −− → EM V and − −− → EM S represent the endmember vectors of vegetation and bare soil/litter in a 2-dimensional space, respectively. The higher value of ESM 2 indicates a better endmember combination for spectral unmixing. It is obvious that the value of ESM 2 is proportional to the area of the triangle constructed by the three endmembers. The best endmember combination is, in effect, corresponding to the largest area of the triangle among all the triangles constructed by any three vertices. Considering that the water endmember is fixed to the pixel nearest the coordinate origin, we just need to search for the other two endmembers (vegetation and soil/litter endmembers) from the remaining pixels so that the three vertices construct the largest triangle. To accelerate the process of endmember identification, the convex hull Graham's scan algorithm [41] that was demonstrated to work well in Yang et al. [42] was applied. The assumption behind this process is that the three pixels that construct the largest triangle area must be located at the convex hull of all the pixels scattered in a two-dimension space. The convex hull was first extracted for the spectral scatter plot by implementing Graham's scan [41], which is a time-efficient method of computing the convex hull of a finite set of points in the plane.
Since the vertex closest to the coordinate origin is not exactly the origin itself (Figure 4b), the origin cannot represent the endmember of water. We thus initiated coordinate translation to ensure that the vertex closest to the coordinate origin matches the origin before the implementation of linear spectral unmixing ( Figure 5). The coordinate translation was done in such a way that it would not impact the triangle shape made by the endmembers; therefore, it would not affect the unmixed fractions. As a result, the water endmember was moved to the coordinate origin in the spectral mixing space. The vertex of the pixel cloud triangle close to the axis of the NIR band represents the vegetation endmember, while the one close to the axis of the red band represents the bare soil/litter endmember.

Linear Spectral Unmixing Model
As mentioned above, linear spectral unmixing was applied to obtain the fractions of vegetation, soil/litter, and water for any pixel in an image. Linear spectral unmixing in the translated spectral mixing space is the process of solving the following linear equations: where ρ (RED) and ρ (NIR) are the spectra of the mixed pixel in the red and NIR band; f v , f s , and f w are the fractions of vegetation, bare soil/litter, and water, respectively; ρ V(RED) and ρ V(NIR) denote the spectra of vegetation endmember in the red and NIR band, whereas ρ S(RED) and ρ S(NIR) represent the spectra of bare soil/litter endmember. Since the coordinate origin is the water endmember, its spectra (ρ W(RED) , ρ W(NIR) ) were assigned as (0, 0). Equation (2) can be rewritten as: In a typical grassland scene, the sum of vegetation, bare soil/litter, and water fractions should be equal to one. The fraction of water can thus be calculated as: Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 20 fractions. As a result, the water endmember was moved to the coordinate origin in the spectral mixing space. The vertex of the pixel cloud triangle close to the axis of the NIR band represents the vegetation endmember, while the one close to the axis of the red band represents the bare soil/litter endmember.

Linear Spectral Unmixing Model
As mentioned above, linear spectral unmixing was applied to obtain the fractions of vegetation, soil/litter, and water for any pixel in an image. Linear spectral unmixing in the translated spectral mixing space is the process of solving the following linear equations: where ρ(RED) and ρ(NIR) are the spectra of the mixed pixel in the red and NIR band; fv, fs, and fw are the fractions of vegetation, bare soil/litter, and water, respectively; ρV(RED) and ρV(NIR) denote the spectra of vegetation endmember in the red and NIR band, whereas ρS((RED) and ρS(NIR) represent the spectra of bare soil/litter endmember. Since the coordinate origin is the water endmember, its spectra (ρW(RED), ρW(NIR)) were assigned as (0, 0). Equation (2) can be rewritten as: In a typical grassland scene, the sum of vegetation, bare soil/litter, and water fractions should be equal to one. The fraction of water can thus be calculated as:

Refinement
In addition to Equation (4), a valid spectral unmixing requires that each endmember fraction is

Refinement
In addition to Equation (4), a valid spectral unmixing requires that each endmember fraction is within the range of zero to one. The translated pixel cloud triangle in the 2-dimensional spectral mixing space is shown in Figure 6. EM V and EM S are the pair of pixels with the highest ESM 2 value. It is apparent that the pixels within the triangle with the vertices of EM V , EM S , and EM W have a meaningful fraction of vegetation, bare soil, and water (i.e., from 0 to 1). However, not all the pixels are always located within this triangle (Figure 4b). For any outlier, the proposed approach will search for the point on the boundary of this triangle with the minimum Euclidean distance to the outlier. In particular, two zones were used to deal with the outlier pixels ( Figure 5). For the pixels in Zone 1, the unmixing solution is the closest endmember (i.e., EM V or EM S ), indicating they are pure pixels of vegetation or bare soil class. On the other hand, the solution for the pixels in Zone 2 is the closest point on either side of the triangle. From this point of view, the Refinement was able to ensure the validity of unmixed fractions for all the pixels in the entire image for each year. data we downloaded from the nearby weather station confirm that the two years had distinctly different precipitation conditions. The early summer of 2005 was very dry, with total precipitation of 36.8 mm in May, 96.2 mm in June, and 16.8 mm from 1-14 July (i.e., right before the image was acquired). On the contrary, the early summer of 2011 was wetter, and its total precipitation was 93.2 mm in May, 77.3 mm in June, and 61.7mm 1-14 July. Within the three days before the images were taken, the study area had 1.8 mm of rainfall on July 13 of 2005, but with 52.2 mm on 13 July 2011.  Figure 7. Visual interpretation clearly indicates low green vegetation cover (purple in Figure 7a), high bare soil/litter cover (dark brown in Figure 7b), and low water content (light blue in Figure 7c) in badlands. On the contrary, there is obvious high green vegetation cover (green in Figure 7a) and low bare soil/litter cover (light yellow in Figure 7b) along the Frenchman River running from northwest to southeast. These results suggest that the unmixing approach resulted in a reasonable FGVC estimation for the study area. However, water content is notably low (light blue to white in Figure  7c) along the Frenchman River, which is likely due to high green vegetation cover in the form of the dense green shrubs and tall grasses growing along the riverbank.

Vegetation Fraction Estimation and Accuracy Assessment
Linear regression between the endmember fraction and the field observed green vegetation cover is then implemented to evaluate the spectral unmixing results. The field data are randomly split into 70% and 30% groups, with 70% to form the regression between unmixed FGVC and observed green cover and 30% for further calibration and evaluation. Linear regression is used because many previous studies have demonstrated its usefulness for evaluating spectral unmixing through its parameters, including the coefficient of determination (R 2 ), mean error (ME), mean absolute error (MSE), and root mean square error (RMSE). A high R 2 , as well as low MAE and RMESE, indicates a better quality of spectral unmixing. After validation, the vegetation maps over multiple years were produced, and vegetation change maps and pixels with significant changes (p < 0.01 or p < 0.05) between 1999 to 2014 were calculated. ME = 1 n Σ(X predicted − X observed ), where X predicted and X observed are the predicted or observed vegetation fraction cover, and n is the number of samples used to validate the predicted vegetation cover.

Time Series Analysis
Following the time series analysis conducted by Bonney et al. [43], we used the unmixed FGVC maps for trend analysis using the "raster" package in R [44]. The changes in the FGVC magnitude were mapped on a per pixel basis through linear regressions, where FGVC was the dependent variable and the image acquisition year was the independent variable. Specifically, the slope of each linear regression was mapped for each pixel to indicate the FGVC change magnitude. The significance of per pixel change was also mapped using the corresponding p-values of a t-test where p values were mapped as p < 0.01, p = 0.01 to 0.05, and p > 0.05. FGVC trends with p ≤ 0.05 were considered a significant change in green vegetation, and FGVC trends with p > 0.05 were considered an insignificant change over the years.
The compiled FGVC maps were also used to extract FGVC data for each of the three dominant landscape units (i.e., badlands, uplands, and riparian zone) in the study area. The effects of precipitation from last 3 days, 1 week, 2 weeks, 3 weeks, and 1 month on vegetation were examined in these three landscape units through a Pearson's r correlation analysis. The significant r values suggested a pronounced impact of the precipitation on the green vegetation cover over time.

Unmixing Results
A total of 15 Red-NIR spectral feature spaces (i.e., scatter plots) were constructed for the 15 images under study, and the two example feature spaces with one from 14 July 2005, and the other from 15 July 2011, are displayed in Figure 6. The overall form of the image-bases spectral spaces is consistent between the two years, although the distributions of the mixed pixel values within the feature space triangles vary ( Figure 6). The similarity of the spectral spaces, as well as those not shown in the figure Figure 7. Visual interpretation clearly indicates low green vegetation cover (purple in Figure 7a), high bare soil/litter cover (dark brown in Figure 7b), and low water content (light blue in Figure 7c) in badlands. On the contrary, there is obvious high green vegetation cover (green in Figure 7a) and low bare soil/litter cover (light yellow in Figure 7b) along the Frenchman River running from northwest to southeast. These results suggest that the unmixing approach resulted in a reasonable FGVC estimation for the study area. However, water content is notably low (light blue to white in Figure 7c) along the Frenchman River, which is likely due to high green vegetation cover in the form of the dense green shrubs and tall grasses growing along the riverbank. A quantitative evaluation was then conducted to validate spectral unmixing results using ground observation data. A strong linear relationship (R 2 = 0.86) between the unmixed FGVC and 70% of observed green cover data is observed (Figure 8a). The regression-derived relationship between unmixed green vegetation cover and 70% of observed field data was applied to the entire A quantitative evaluation was then conducted to validate spectral unmixing results using ground observation data. A strong linear relationship (R 2 = 0.86) between the unmixed FGVC and 70% of observed green cover data is observed (Figure 8a). The regression-derived relationship between unmixed green vegetation cover and 70% of observed field data was applied to the entire map of the unmixed FGVC, which was then validated by the remaining 30% of observed green vegetation cover data. The results indicate that the correction led to a strong FGVC estimation with an R 2 value of 0.93 and an RMSE of 7.05% (Figure 8b). A quantitative evaluation was then conducted to validate spectral unmixing results using ground observation data. A strong linear relationship (R 2 = 0.86) between the unmixed FGVC and 70% of observed green cover data is observed (Figure 8a). The regression-derived relationship between unmixed green vegetation cover and 70% of observed field data was applied to the entire map of the unmixed FGVC, which was then validated by the remaining 30% of observed green vegetation cover data. The results indicate that the correction led to a strong FGVC estimation with an R 2 value of 0.93 and an RMSE of 7.05% (Figure 8b).

Time Series Results
The FGVC maps from 1999 to 2014 are displayed in Figure 9. The field data that were reserved for validation were also used to calculate the RMSE for each year between the unmixing FGVC and observed green vegetation cover. The RMSE ranges from the lowest value of 3.35% in 2010 to the highest 22.33% in 2003. It is clear that the farther the image acquisition date was away from the field data collection season (15 June to 15 July), the greater the RMSE is. For example, the unmixing result in 2003 has a large error, while the 2013 image was taken on August 10, about a month after the field data were collected. In contrast, the FGVC in 2010 has a small RMSE, corresponding to a complete match between the field season and image acquisition date of 26 June.
The amount of green vegetation varies substantially from year to year, and the change is not uniform across space (Figure 9). To further investigate the difference in the greenness in the study area, we selected the unmixed FGVC maps from 14 July 2005 and 15 July 2011, for a detailed comparison. Consistent with higher precipitation in 2011 than in 2005, the unmixed FGVC has a generally higher vegetation cover in the year of 2011 (i.e., the wet year) than in the year of 2005 (i.e., the dry year). This is demonstrated by more areas with a light green color (i.e., closer to 100% FGVC) and by badlands having a lighter purple color (i.e., farther from 0% vegetation) in 2011 in comparison

Time Series Results
The FGVC maps from 1999 to 2014 are displayed in Figure 9. The field data that were reserved for validation were also used to calculate the RMSE for each year between the unmixing FGVC and observed green vegetation cover. The RMSE ranges from the lowest value of 3.35% in 2010 to the highest 22.33% in 2003. It is clear that the farther the image acquisition date was away from the field data collection season (15 June to 15 July), the greater the RMSE is. For example, the unmixing result in 2003 has a large error, while the 2013 image was taken on August 10, about a month after the field data were collected. In contrast, the FGVC in 2010 has a small RMSE, corresponding to a complete match between the field season and image acquisition date of 26 June.
The amount of green vegetation varies substantially from year to year, and the change is not uniform across space (Figure 9). To further investigate the difference in the greenness in the study area, we selected the unmixed FGVC maps from 14 July 2005 and 15 July 2011, for a detailed comparison. Consistent with higher precipitation in 2011 than in 2005, the unmixed FGVC has a generally higher vegetation cover in the year of 2011 (i.e., the wet year) than in the year of 2005 (i.e., the dry year). This is demonstrated by more areas with a light green color (i.e., closer to 100% FGVC) and by badlands having a lighter purple color (i.e., farther from 0% vegetation) in 2011 in comparison with 2005 images (Figure 10a,b). The enlarged maps (Figure 10(a1,b1) further show that a few ponds full of water in the 2011 FGVC map (purple colored polygons in Figure 10(b1)) were occupied by dense green vegetation (dark green colored polygons in Figure 10(a1)) in the 2005 FGVC map. These results confirm that water availability is the primary driver of variations in green vegetation cover in this semi-arid area.
with 2005 images (Figure 10a,b). The enlarged maps (Figure 10(a1,b1) further show that a few ponds full of water in the 2011 FGVC map (purple colored polygons in Figure 10(b1)) were occupied by dense green vegetation (dark green colored polygons in Figure 10(a1)) in the 2005 FGVC map. These results confirm that water availability is the primary driver of variations in green vegetation cover in this semi-arid area.  FGVC trends from the year 1999 to 2014 were mapped over the entire study area (Figure 11a). An increase in FGVC of less than 1% per year is observed in most of the study area. Only a small portion of the study area has experienced a decreasing trend, mainly along the Frenchman River and its tributaries, with a reduced rate of less than 1% per year. Significant changes, both increasing and FGVC trends from the year 1999 to 2014 were mapped over the entire study area (Figure 11a). An increase in FGVC of less than 1% per year is observed in most of the study area. Only a small portion of the study area has experienced a decreasing trend, mainly along the Frenchman River and its tributaries, with a reduced rate of less than 1% per year. Significant changes, both increasing and decreasing (Figure 11b), occurred mainly in the west half of the study area or along the river. FGVC trends from the year 1999 to 2014 were mapped over the entire study area (Figure 11a). An increase in FGVC of less than 1% per year is observed in most of the study area. Only a small portion of the study area has experienced a decreasing trend, mainly along the Frenchman River and its tributaries, with a reduced rate of less than 1% per year. Significant changes, both increasing and decreasing (Figure 11b), occurred mainly in the west half of the study area or along the river. The unmixed FGVC for badlands, uplands, and riparian zone were further extracted from each FGVC map, and their temporal variations are displayed in Figure 12. FGVC is always higher in the riparian zone, followed by the uplands and then the badlands. Consistent with the significance map (Figure 11b), the linear regression trend lines in Figure 12a also suggest a slight increase in FGVC in the uplands and badlands but a decrease in the riparian zone. The vegetation greenness in badlands and uplands are most highly correlated with the total precipitation accumulated from 3 days before the images were taken (Table 2). Similarly, the last 3-day total precipitation also showed a slightly increasing trend during the study period (Figure 12b). On the other hand, vegetation greenness in The unmixed FGVC for badlands, uplands, and riparian zone were further extracted from each FGVC map, and their temporal variations are displayed in Figure 12. FGVC is always higher in the riparian zone, followed by the uplands and then the badlands. Consistent with the significance map (Figure 11b), the linear regression trend lines in Figure 12a also suggest a slight increase in FGVC in the uplands and badlands but a decrease in the riparian zone. The vegetation greenness in badlands and uplands are most highly correlated with the total precipitation accumulated from 3 days before the images were taken (Table 2). Similarly, the last 3-day total precipitation also showed a slightly increasing trend during the study period (Figure 12b). On the other hand, vegetation greenness in the riparian zone is less impacted by the short-term precipitation, as shown by their weaker correlation in comparison with the correlations for the badlands and uplands. Vegetation in the riparian zone tends to have a stronger correlation with the last 7-day precipitation ( Table 2). the riparian zone is less impacted by the short-term precipitation, as shown by their weaker correlation in comparison with the correlations for the badlands and uplands. Vegetation in the riparian zone tends to have a stronger correlation with the last 7-day precipitation ( Table 2).
. Figure 12. The unmixed FGVC extracted from typical badlands, typical riparian zone, and typical uplands from 1999 to 2014, their linear regression trend lines (a), and the total precipitation from the last three days to the last month before the images were taken and the last three-day total precipitation linear regression trend line (b).

Discussion
As an important branch of soft classification, spectral unmixing has been commonly used in many remote sensing applications, especially for medium/low-resolution imagery. In recent years, several studies have started to unmix fine-scale imagery (i.e., IKONOS and WorldView-2) through the analysis of spectral mixing space [36,37]. In comparison with the conventional spectral unmixing, these studies identified and determined the endmembers by investigating the shape of pixel clouds in spectral mixing spaces. Specifically, the shape of the pixel cloud in the spectral mixing space constructed by the red and NIR band was found to be a triangle, and the three vertices of this triangle were 100% vegetation, bare soil/litter, and water endmembers. The image-based red-NIR space feature space is consistent over the years, and similarly, vegetation endmembers are also stable for this study area, suggesting a simple linear mixture is able to construct a consistent grassland reflectance. Using the endmembers identified from each image, spectral unmixing was successfully conducted, and the unmixed FVGC with an acceptable RMSE indicates fully automated spectral Figure 12. The unmixed FGVC extracted from typical badlands, typical riparian zone, and typical uplands from 1999 to 2014, their linear regression trend lines (a), and the total precipitation from the last three days to the last month before the images were taken and the last three-day total precipitation linear regression trend line (b).

Discussion
As an important branch of soft classification, spectral unmixing has been commonly used in many remote sensing applications, especially for medium/low-resolution imagery. In recent years, several studies have started to unmix fine-scale imagery (i.e., IKONOS and WorldView-2) through the analysis of spectral mixing space [36,37]. In comparison with the conventional spectral unmixing, these studies identified and determined the endmembers by investigating the shape of pixel clouds in spectral mixing spaces. Specifically, the shape of the pixel cloud in the spectral mixing space constructed by the red and NIR band was found to be a triangle, and the three vertices of this triangle were 100% vegetation, bare soil/litter, and water endmembers. The image-based red-NIR space feature space is consistent over the years, and similarly, vegetation endmembers are also stable for this study area, suggesting a simple linear mixture is able to construct a consistent grassland reflectance. Using the endmembers identified from each image, spectral unmixing was successfully conducted, and the unmixed FVGC with an acceptable RMSE indicates fully automated spectral unmixing is possible for the semi-arid mixed grasslands and yields reasonable vegetation fraction estimation.
Caution should be taken when applying the automated spectral unmixing approach, because a close look at the unmixing images suggests a systematic FGVC underestimation (Figure 8a), resulting from the overestimation of water percentage in the mixed pixels. A potential explanation for water overestimation is the contribution of standing vegetation shadows to the unmixing results but not to the field observation data. In the field, the shadows of standing grasses would not contribute to the visual interpretation of the green vegetation cover, and thus the canopy cover will be estimated properly by a field crew. However, the shadows of grasses would contribute to spectral unmixing because shadow spectra are inherent in each mixed pixel, and the shadow spectra are similar to water given lower spectral values of shadow in the red and NIR region. Similarly, Li et al. [45] and Jiang et al. [46] also found that the presence of shadow leads to NDVI saturation, which then results in the overestimation of FVGC.
It is also worth noting that this study only focused on green vegetation cover, so it was not an issue to combine bare soil and litter as one endmember. However, if the goal is to obtain litter abundance or soil coverage, SWIR or the decay pigment index may be added to separate litter and soil [30,47]. Still, attention should be paid to ensure that these endmembers are not correlated. Otherwise, the inversion of the linear spectral mixture model would become unstable, and the estimated fractions would be highly sensitive to random error [48,49]. It has been advised that an appropriate spectral mixing space is required to avoid endmember collinearity and multi-collinearity, to further improve the quality of spectral unmixing [37].
The study showed an increasing trend in FGVC in general for badlands and uplands in the study area. The results are consistent with previous studies; for example, Ju and Masek [50] found a greening trend from 1984 to 2012 using Landsat in Canadian grasslands in Southern Saskatchewan. In our study site, however, the FGVC in the riparian zone has shown a slight negative trend over the years. The different trends in different landscape units are likely a result of different plant communities and their various sensitivities to climatic conditions, especially precipitation. Native grass-dominated uplands and badlands are more likely dependent on precipitation (Table 2), while the riparian zone with more tall shrubs and dense invasive vegetation is not as sensitive to the precipitation as grass, especially to 3-day precipitation. Similarly, Flanagan etc. also suggested that grasslands respond greatly to annual precipitation variation [51], and they further indicated that the responses are not even. Knapp and Smith found that grasslands experience an increase in productivity in wet years and a reduction during drought years [52].
Only a small portion of the pixels exhibited a significant change in our study area. This is understandable because the area represents the Prairie Grasslands natural region in Canada and has been protected from human activities for decades. The pixels undergoing significant change are likely a result of climate dynamics and natural disturbances. Other than precipitation, moisture has been identified as a limiting factor for plant growth in semi-arid ecosystems like the northern mixed-grass and shortgrass prairie [8], and litter is a critical structural component which conserves moisture. Two other potential factors may have contributed to a decreased green vegetation cover in the riparian zone by reducing litter amount and further leading to lower soil moisture. The first is that the GNP reintroduced bison into the park in 2005 [12], and the second is a massive wildfire that occurred in the riparian zone of the study area in 2013 [8]. However, these are only possible explanations, and future work is warranted to confirm this. In specific, longer-term changes (e.g., the past 40 years) and the reasons behind these significant changes need to be further explored.
Other than the different drivers that may contribute to the FGVC trends, the phenological stage when an image is taken could also impact the FGVC values in certain years and further contribute to the trends. We have five images taken from August (1999,2003,2008,2012,2014), and apparent lower FGVC values in 1999, 2008, and 2012 ( Figure 12) could have been a combined result of climate conditions and plants experiencing a different phenological stage. Sonnenschein et al. (2011) also found high variability in vegetation temporal trends in grasslands, and they further indicated that high phenological variability potentially affects a single year's relationship between different vegetation estimates and leads to varying trends. They suggested that stratifying a study region prior to trend analyses may address the varying trends. Following the advice, we performed a trend analysis over three different landscape units in the study area and did observe a different trend between units. Still, when interpreting our trend analysis results, the readers should keep in mind that phenological stages of varying vegetation communities may have played a role in our trend analysis results.

Conclusions
In this study, an automated spectral mixture approach is applied to a series of Landsat images to unmix green vegetation cover for a mixed grassland between 1999 and 2014. The spectral unmixing results correctly reflected the spatial and temporal variations of green vegetation cover in the study area. An increasing trend in FGVC is observed in general for the study area during the period, and this increase is mainly observed in the badlands and uplands, as a result of temporal variation in the short term precipitation. In contrast, a decreasing trend is observed in riparian zone vegetation, likely due to its insensitivity to precipitation. Further research should explore the ways to remove the impact of grass shadows in the spectral unmixing, distinguish bare soil and litter by introducing another endmember, and understand the drivers behind the pixels with a significant decrease in this endangered grassland environment. With precipitation in Canadian prairies predicted to fluctuate under a warmer climate [53], further studies on this topic will be essential to understand the implications of climate dynamics on grassland functions [52].