Optimizing Spectral Libraries from Landsat Imagery for the Analysis of Habitat Richness Using MESMA

Abstract

Spectral mixture analysis of satellite images, such as MESMA (multiple endmember spectral mixtures analysis), can be used to obtain fraction images in which the abundance of each land occupation class is represented at the pixel level, which is crucial for the analysis of heterogeneous landscapes in which types of habitats vary at fine spatial scales. The objective of this work is to analyze the influence of spectral libraries of various characteristics on the performance of MESMA. To this end, eight spectral libraries from Landsat satellite images were elaborated with different characteristics in terms of size, composition, and temporality. The spectral libraries were optimized using the iterative selection of endmembers (IES) method with the MESMA technique to obtain the fraction images considering five habitat classes (forest, shrubland, grassland, water, and rock and bare soil). The application of MESMA resulted in the classification of more than 95% of pixels in all cases with a root mean square error (RMSE) less than or equal to 0.025. Validation of the fraction images through linear regressions resulted in an RMSE ≥ 0.35 for the shrubland and grassland classes, with a lower RMSE for the remaining classes. A significant influence of library size was observed, as well as a significant effect of temporality, with the best results obtained for the largest monotemporal libraries.

1. Introduction

Remote sensing allows for the identification of Earth’s surface characteristics by obtaining information from the electromagnetic radiation reflected from the various materials on the terrestrial surface [1]. The need to understand ecosystems on a global scale and in a broad temporal context has resulted in the inclusion of remote sensing in numerous ecological studies on various topics involving the use of satellite images for the analysis of environmental characteristics [2].

The spatial, radiometric, spectral, and temporal scale at which the data derived from satellite imagery are available varies, depending on the type of sensor from which they are obtained (airborne or satellite; active or passive; multispectral or hyperspectral). The consideration of spatial, radiometric, spectral, and temporal characteristics is an important step in the application remote sensing to the study of landscapes [3]. On the other hand, the availability of remote sensing data determines, in many cases, the study approaches, as the complexity of the landscape and the working scale determine the optimal type of data and methodological approaches necessary to conduct studies [4]. Multispectral images from sensors, such as MODIS, Landsat and Sentinel-2 data, have been widely used to conduct studies at medium–fine spatial scales [5,6]; however, when the approach involves temporal analysis of the landscape, the Landsat sensor, the data of which have been obtained since the early 1970s, is the most widely used [7].

In recent decades, several methods enabling the classification of multispectral satellite images have been developed. Among them, we can differentiate between (I) hard classifications that assign a single class to each pixel and are commonly implemented through the application of maximum likelihood algorithms, neural networks, support vector machines (SVMs), or decision trees and (II) soft classifications, such as those resulting from the application of spectral mixture analysis (SMA), which takes into account the existence of spectral mixtures within a pixel, providing results on the abundance of each class within a single pixel [3,8]. The selection of one of these classification strategy affects the accuracy of the resulting spatial models. Therefore, the spatial complexity of the target landscape determines the optimal scale at which the analysis performs and therefore the most appropriate method to carry out image classification [3]. It has been observed that soft classification methods allowing for the analysis of class abundances at the pixel level have achieved better results in fragmented and heterogeneous landscapes than hard classification, as they consider the variability of existing classes at this level [9]. These types of complex landscapes present a high number of fine-scale habitats that classic image classification methods cannot model. Moreover, with to analysis of their diversity at pixel scales, spectral mixture methods (also known as spectral unmixing) have achieved satisfactory results in obtaining measures of diversity, although their capacity to accurately model numerous landcover classes remains limited [10].

Spectral unmixing techniques are based on the use of spectral signatures that represent the different classes in a study area. These spectral signatures can be obtained in the laboratory from reference materials [11], measured in the field [12], through radiative transfer models [13], or from the satellite images used for classification [14]. Then, the spectral signatures are included in a library, which is used to perform the spectral unmixing of each pixel, obtaining the abundance data for each existing class within a pixel.

As observed in previous studies [15], the quality of the results obtained through this type of method depends directly on how well the spectral signatures included in the libraries represent the existing classes in the studied scene. Therefore, among the most important factors when applying this classification method is the use of an optimal spectral library in terms of size and composition. To improve the quality of the libraries, spectral signature selection techniques have been developed to include in the libraries only the optimal spectral signatures that best represent the existing classes in the study area [16]. Among the most used techniques to optimize spectral libraries are EAR (endmember average root mean square error), which selects spectral signatures that provide the lowest RMSE after performing a linear SMA between the signatures of a class [17]; MASA (minimum average spectral angle), which selects the spectral signatures that present the minimum value for the spectral angle [18]; CoBI (count-based spectral signature selection index), which selects the spectral signatures that model a greater number of spectra within each class [19]; and IES (iterative selection of spectral signatures), which selects the group of spectral signatures that produces the highest value of the Kappa coefficient [16]. Among these techniques, IES provides the best results, as, in addition to considering the variability of spectral signatures within a given class, it also considers the differences between classes [20].

The linear SMA technique allows for the spectral unmixing of pixels with a reduced number of endmembers and does not consider the temporal and spatial variability between the spectral signatures themselves [9]. However, in recent years, more advanced spectral mixture analysis techniques have been developed, such as MESMA (multiple endmember spectral mixture analysis), allows for the analysis to be carried out with a greater number of spectral signatures and analysis of each pixel with a smaller independent group of spectral signatures [21]. In MESMA, unlike traditional spectral unmixing approaches, the spectral variability within each class is considered, and the pixels are modelled using a set of spectral signatures that is not fixed and varies throughout the image. In this case, the accuracy of the results can be improved and compared to traditional approaches in which such variation is not considered, and a fixed set of spectral signatures is used to model the full image, resulting in possible inaccurate estimates of coverage fractions [22,23]. MESMA allows for spectral unmixing using a combination of two or more spectral signatures on each modelled pixel, including a shade fraction caused by the topography, in order to assume the effect that it has on the spectral signatures. Finally, it assigns the best-fitting model to each pixel based on the lowest value of RMSE [9]. To date, this technique has been used to model a limited number of occupation classes, usually two or three [13,24,25,26]. However, modelling a larger number of classes in complex landscapes remains a challenge that would require, in the first instance, the use of highly optimized libraries.

The general objective of this work is to analyze the influence of different types of spectral libraries (according to their size, composition, and temporality) on the operation of MESMA to carry out spatiotemporal analysis in a complex landscape characterized by a high number of classes of habitats at a fine spatial scale.

2. Materials and Methods

2.1. Study Area

The study area is located in the northwest of the Iberian Peninsula, coinciding with a rectangle of 110 × 33 km that includes the central–western sector of the Cantabrian Mountains (Figure 1).

The study area comprises the Eurosiberian biogeographic region in the northern part and the Mediterranean region in the southern portion. Mean annual precipitation ranges between 700 and 2200 mm, and mean annual temperature ranges between 4 and 14 °C [27]. The lithology is heterogeneous, the topography is abrupt, and the elevation ranges from 400 to 2400 m in short distances.

The wide range of environmental conditions in the study area favors the existence of heterogeneous landscapes with a high richness of habitats (i.e., a number of different habitats) on a fine spatial scale. The main tree species in the study area belong to genera such as Fagus (F. sylvatica L.) and Quercus (Q. robur L., Q. petrea Matt., Q. ilex L., and Q. pyrenaica Willd.) in high-altitude areas, as well as birch trees, deciduous forests, and pine plantations, with species such as Pinus sylvestris L., P. nigra Arn., and P. radiata D. Don. In terms of shrub vegetation, heaths and bushes are prominent, dominated by species of families such as Ericaceae, Fabaceae, and Cistaceae. Pastures and meadows that intermingle with wooded and shrub areas, as well as peaks dominated by rocky surfaces and natural grasslands. The study area also contains water bodies, such as the Luna, Porma, and Riaño reservoirs in the province of Leon.

According to the classes of habitats identified in the study area in previous work [10], five classes were considered for elaboration in this study: shrubland (S), forest (F), grassland (G), water (W), and rock and bare soil (R). These types of habitats were visually identified in available orthophotographs of the study area (Figure 2). The habitats intermingle at a fine scale across the entire study area, with a predominance of forests on the northern face of the mountain range.

2.2. Remotely Sensed and Ancillary Data

Four Landsat satellite images from 1990, 2000, 2010, and 2020 were used. The images belong to Landsat collection 2 Level-2 scene with bottom of atmosphere (BOA) reflectance values downloaded from the United States Geological Survey server (USGS: https://earthexplorer.usgs.gov/ accessed on 14 September 2022). The first three images (1990, 2000, and 2010) were obtained from the Landsat 5 TM sensor and were atmospherically corrected with the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm, whereas the 2020 image was obtained from the Landsat 8 OLI sensor and was atmospherically corrected with the Land Surface Reflectance Code (LaSRC) algorithm.

Landsat images were at 30 m pixel in size and georeferenced to the WGS 84/UTM Coordinate Reference System (CRS) 30 N zone (EPSG Code: 32630). All images were cloud-free, and the acquisition dates corresponded to the end of the summer period, when the snow cover is testimonial. Data were obtained from the multispectral bands corresponding to blue, green, red, near-infrared (NIR), and short-wave infrared 1 and 2 (SWIR 1 and SWIR 2). Moreover, we used the closest available orthophotographs to the dates of Landsat imagery at ≤ 50 cm spatial resolution as reference data. Specifically, orthophotographs of the years 1973–1986, 1997–2003, 2008–2011, and 2020 were obtained from the National Aerial Orthophotography Project (PNOA) (IGN; https://www.ign.es/web/ign/portal, accessed on 1 January 2022) and from the Agrarian Technological Institute of Castile and Leon (ITACYL; http://orto.wms.itacyl.es/WMS, accessed on 1 January 2022). All the orthophotographs used in the present study completely covered the study area.

2.3. Methods

The methodology comprises a total of 4 steps: (1) image preprocessing, (2) creation and optimization of spectral libraries, (3) MESMA procedure, and (4) validation (Figure 3). QGIS software [28] was used for image processing, spectral library creation and optimization, the MESMA procedure, and validation. Statistical analyses were performed with RStudio [29].

2.3.1. Image Preprocessing

Landsat images were converted to reflectance using the scale and offset factors recommended by the USGS. Subsequently, anomalous values in the various bands were eliminated, assigning a value of 1 to pixels that presented values greater than 1. The six bands with which the analysis would be carried out for each of the satellite images were stacked, resulting in four groups of bands (one per year under study: 1990, 2000, 2010, and 2020). These groups of bands are referred to as “satellite images” throughout this work. Starting satellite images were projected to the ETRS89 UTM zone 30N coordinate reference system. Subsequently, each of the groups of bands was cut out from a 110 × 33 km rectangular polygon defining the study area to obtain the satellite images for subsequent analysis.

2.3.2. Creation and Optimization of Spectral Libraries

Spectral libraries were elaborated according to the following criteria: (I) size of the library (large or small), (II) composition according to the type of sensor from which the spectral signatures were extracted (Landsat TM or Landsat OLI), and (III) temporality (monotemporal or multitemporal). In the multitemporal libraries, spectral signatures of the various sensors were combined from the images corresponding to the four years under study, and the monotemporal libraries were made from spectral signatures of a single year and therefore of the same sensor.

Prior to MESMA analysis, it is necessary to define the endmembers (pure spectral signatures) that represent each class of habitat defined in the study area and that will form part of the spectral libraries used for spectral unmixing of the pixels of the satellite images [21]. In our study, the endmembers were extracted from the used satellite images. Although their spectral signatures were extracted from Landsat imagery, these pure areas were identified through the reference orthophotographs as areas 90% covered by a single class. Then, we implemented a stratified random design to distribute 250 30 × 30 m plots (Figure 4). Within each habitat class, distinguishing between forest, shrubland, grassland, rock and bare soil, and water, a total of 50 plots were distributed in areas that did not change over the study period according to the orthophotographs. To make up for the lack of coverage of some orthophotographs, combinations of Landsat bands were visualized to identify the types of habitats in the study area as false-color compositions for the identification of healthy and well-developed vegetation in front of water areas and bare soil. This combination of bands was used exclusively to select a plot corresponding to a pure zone.

From the sample of 250 30 × 30 m plots, libraries were generated with the “Spectral Library Tool” plugin for QGIS 3.14 software [30]. The input data for the elaboration of these spectral libraries were the 250 previously generated 30 × 30 m plots of and the satellite image for each year. The operation of the indicated plugin consists of generating the centroid of the inserted polygons and extracting the reflectance value of each pixel for each of the bands of the satellite images. The generated spectral libraries are described below.

Two multitemporal libraries with inequitable distribution between sensors (equitable among years) were made up of 25% spectral signatures of each image, i.e., 75% spectral signatures corresponding to the Landsat 5TM sensor obtained from the images for the years 1990, 2000, and 2010 and 25% spectral signatures corresponding to the Landsat 8OLI sensor obtained from the image for the year 2020. The first library, Library 1 (L1), was the largest and was generated using all plots and all years (250 spectral signatures per year; 1000 spectral signatures), and the second library, Library 3 (L3), was smaller and was generated using only 100 of the 250 plots per year, with 20 plots randomly selected per habitat type (100 spectral signatures per year; 400 spectral signatures).

Another two multitemporal libraries with equitable distribution between sensors were made up of 50% spectral signatures of the Landsat 5TM sensor (years 1990, 2000, and 2010) and 50% spectral signatures of the Landsat 8OLI sensor (year 2020). The first library, Library 2 (L2), was the largest and was generated using 250 plots of the Landsat 8OLI sensor and by reducing the 750 remaining spectral signatures of the Landsat 5TM sensor to 250 plots, and the second library, Library 4 (L4), was smaller and was generated using 100 of the 250 plots of the Landsat 8OLI sensor and by reducing the 300 remaining spectral signatures of the Landsat 5TM sensor to 100 plots.

The monotemporal libraries were elaborated from the 250 spectral signatures of each year extracted for the elaboration of L1. A total of four monotemporal libraries were obtained (for the years 1990, 2000, 2010, and 2020), with 250 spectral signatures each. Ultimately, a total of eight spectral libraries of varying sizes and compositions were generated, differentiating between multitemporal libraries (L1, L2, L3, and L4) and monotemporal libraries (1990, 2000, 2010, and 2020).

The spectral libraries were optimized to extract the most representative groups of spectral signatures of each class from each library for subsequent image classification. This step reduces the cost of executing MESMA and improves the accuracy of the results [16]. To this end, the iterative selection of endmembers (IES) technique was used, which has yielded satisfactory results in previous studies in terms of improving the accuracy of image classification [16,31]. IES is an automatic method that extracts the group of spectral signatures that produces the highest value of Cohen’s Kappa coefficient [32] by iteratively adding and removing spectral signatures from the analysis. When an endmember is included that does not improve the coefficient, no more endmembers are included in the group of signatures. For the application of the IES technique, 0 and 1 were assigned as minimum and maximum allowed fractions, respectively, with a maximum mean square error (RMSE) of 0.025 [16,21,33], obtaining spectral libraries that were smaller than the initial libraries (

Table 1. Characteristics and composition of the developed spectral libraries.

Library Group	Spectral Library	Composition (Sensor)	Composition (Year)	Size (Endmembers)	Size after IES (Endmembers)
MULTITEMPORALWITH INEQUITABLE DISTRIBUTION	L1	75% TM 25% OLI	25% 1990 25% 2000 25% 2010 25% 2020	1000	35
	L3			400	24
MULTITEMPORAL WITH EQUITABLE DISTRIBUTION	L2	50% TM 50% OLI	50% 1990, 2000 & 2010 50% 2020	500	27
	L4			200	14
MONOTEMPORAL	1990	100% TM	100% 1990	250	14
	2000		100% 2000		13
	2010		100% 2010		13
	2020	100% OLI	100% 2020		17

).

2.3.3. MESMA Procedure: Obtaining Fraction Images

The eight optimized spectral libraries and the processed Landsat images were used to perform the MESMA analysis of the four study years (i.e., four satellite images), resulting in 20 MESMA classifications (monotemporal libraries were only applied to their corresponding satellite image). The MESMA plugin [34] was used to perform MESMA analysis, which is based on VIPER Tools 2.0 developed by the VIPER laboratory at UC Santa Barbara.

Pixel modelling was carried out by checking three types of models: one type of habitat + shade fraction, two types of habitats + shade fraction, and three types of habitats + shade fraction, selecting the optimal combination according to the RMSE statistic. The values of the minimum and maximum fraction were set at 0 and 1, respectively, and the maximum shade fraction was restricted to 0.80 [10,25]. A type of photometric shade was included in which the shade has a default reflectance of 0. With respect to the RMSE, a maximum value of 0.025 was assigned so that all pixels in the classification with a value greater than the RMSE would not be modelled.

MESMA analysis resulted in three products: an image of the models in the form of the number of classes of habitats in each pixel, the fraction images for each type of habitat present in a pixel and the shade, and an image with the RMSE value per pixel. The obtained fraction images were normalized by eliminating their shade and maintaining the value of the abundance of each type of habitat for each pixel.

2.3.4. Validation

In order to validate the results obtained from MESMA and to analyze the operation of the generated libraries, we implemented a stratified random sampling design with 50 validation points per year, distributing 10 points in the areas dominated by each class according to the orthophotographs (See Figure A1 in Appendix A).

Then, a 30 × 30 m square buffer was generated around each of the 50 points using QGIS 3.14 software. The 30 × 30 m square polygon was matched to the spatial resolution of the satellite images. Finally, we configurated our reference datasets (one per year) by assigning the fraction of each type of habitat to each 30 × 30 m polygon with the support of the closest orthophotograph. Multiple values of 0.1 were used (Figure 5). On the other hand, we extracted the values resulting from the normalized fraction images resulting from MESMA using the same 50 points described above using the “Point Sampling Tool” of the QGIS software. Then, we performed bivariate linear regressions between each normalized MESMA fraction (predictor variable) and the reference data (response variable). Statistical analysis was performed with Rstudio software, and the RMSE was calculated for each of the regressions performed.

2.3.5. Analysis of the Performance of Spectral Libraries

Mixed linear models were elaborated to determine whether the size and type of library had significant effect on the operation of MESMA, depending on the origin sensor of the spectral signatures and its temporality. The RMSE obtained in the validations for the MESMA normalized fraction images performed with each of the spectral libraries (L1_IES, L2_IES, L3_IES, L4_IES, 1990_IES, 2000_IES, 2010_IES, and 2020_IES) was included as a dependent variable, and the size and temporality/composition of the spectral library were included as independent variables. Finally, the year of the image and the type of habitat were included as random factors.

After elaborating the statistical models, we confirmed that model residuals were normally distributed (using the Shapiro–Wilk test) and homoscedastic (equality of variances graphically displayed). The lmer function of the lme4 package [35] was used to perform linear mixed models, and the plot_model function of the SjPlot package [36] was used to display the results of these models.

3. Results

3.1. Spectral Libraries

The spectral signatures of the optimized spectral libraries generated after applying the IES analysis are shown in Figure 6.

We observed that the optimized spectral libraries included a greater number of endmembers for the classes of shrubland (48 endmembers), grassland (31 endmembers), and rock and bare soil (16 endmembers) than for the classes of water (3 endmembers) and forest (9 endmembers).The reflectance peaks for classes such as forest and shrubland are mainly in the near-infrared range, whereas in the case of the grassland class, peaks appear in the near-infrared and short-wave infrared (SWIR 1) ranges. For the rock and bare soil class, the wavelength in the blue and green bands, i.e., between 0.10 and 0.20 µm, is greater than that of the rest of the classes, which is found between 0 and 0.10 µm. The reflectance peak for this class is in the SWIR 1 range. For the water class, none of the bands has a reflectance value greater than 0.10 µm.

When applying the IES analysis to the various spectral libraries, the Kappa values ranged from 0.17 in all the libraries when a single spectral signature was used to 0.76 with the use of a spectral signature of each class, that is, a total of five spectral signatures. Ultimately, values between 0.90 and 0.95 were obtained among the spectral libraries (

Table 2. Kappa values obtained after IES analysis for each spectral library, the number of spectral signatures included in each library, and the spectral composition of the libraries after optimization with IES.

Library	Kappa 1 Endmember	Kappa 5 Endmembers (1 Per Class)	Final Kappa	Endmembers	Endmembers of Each Class					Sensor	Endmembers %
					W	F	S	G	R
L1_IES	0.17	0.76	0.91	35	-	2	15	3	7	TM	77
					1	-	2	4	1	OLI	23
L2_IES	0.17	0.76	0.91	27	-	1	6	2	1	TM	37
					1	-	6	5	5	OLI	63
L3_IES	0.17	0.78	0.95	24	-	2	10	4	3	TM	80
					1	-	3	1	-	OLI	21
L4_IES	0.17	0.77	0.93	14	-	-	3	-	1	TM	29
					1	1	3	4	1	OLI	71
1990_IES	0.17	0.76	0.93	14	1	1	4	5	3	TM	100
2000_IES	0.17	0.78	0.90	13	1	2	4	3	3	TM	100
2010_IES	0.17	0.88	0.93	13	1	2	3	4	3	TM	100
2020_IES	0.17	0.78	0.95	17	1	2	6	4	4	OLI	100

).

The size of the optimized libraries varied from the 35 spectral signatures in the multitemporal library with the largest unequitable distribution (L1_IES) to 13 spectral signatures in the monotemporal libraries from 2000 and 2010 (2000_IES and 2010_IES, respectively). Furthermore, in the case of the L1_IES, L2_IES, and L3_IES libraries, the difference between the number of spectral signatures of the shrubland class relative to the other classes was greater than in the remaining spectral libraries. In all cases, only one spectral signature of the water class and one or two spectral signatures of the forest class were included.

The multitemporal libraries presented unequal distribution of signatures among sensors. In the case of the L1_IES (77% TM and 23% OLI) and L3_IES (80% TM and 21% OLI) libraries, the percentage of signatures of each sensor was close to the percentage included in the initial L1 and L3 libraries (75% TM and 25% OLI). For the L2_IES and L4_IES libraries, which started as L2 and L4 libraries formed from 50% TM sensor signatures and 50% OLI sensor signatures, respectively, the result was a high percentage of OLI-derived signatures (>60%).

In all cases, the spectral signature corresponding to the water class for the multitemporal libraries was derived from the Landsat OLI sensor.

3.2. MESMA Results and Optimal Spectral Libraries

The size of the optimized spectral libraries used for classification of the satellite images with MESMA influenced the number of models applied to each pixel of the Landsat images; the larger the library, the more models applied. The percentage of pixels classified with an RMSE value of less than 0.025 was greater than 95% in all cases (

Table 3. Percentage of classified pixels after applying MESMA with the different spectral libraries and the percentage of pixels classified by class according to the MESMA analyses performed with the optimized monotemporal libraries (RMSE ≤ 0.025).

		% of Classified Pixels
LIBRARY		1990	2000	2010	2020
L1_IES		98.49	96.95	99.39	97.29
L2_IES		98.77	97.6	99.39	97.71
L3_IES		98.67	97.32	99.01	97.44
L4_IES		97.56	95.33	98.52	95.71
MONOTEMPORAL		97.87	97.78	97.75	97.36
CLASS	Water	0.84	1.47	0.77	0.68
	Forest	21.68	27.85	23.35	33.68
	Shrubland	31.7	24.44	37.48	26.51
	Grassland	29.4	28.09	26.74	26.11
	Rock and bare soil	16.38	18.15	11.66	13.01

). In addition, a an improved performance of the optimized monotemporal libraries was observed.

In the classifications made with the optimized multitemporal libraries, a difference was observed between the results obtained for the images from the years 1990 and 2010, for which the percentage of classified pixels was around 98–99%, and the images from the years 2000 and 2020, for which the percentage of classified pixels was less than 98%. The L2_IES library presented high percentage values of classified pixels for all the images under study. A significant difference was observed in the results obtained for the images for the years 2000 and 2020 after classifications made with the multitemporal library L4_IES, for which the percentage of classified pixels was between 95–96%, which was the lowest percentage of those obtained with the various spectral libraries.

The optimized monotemporal libraries (1990_IES, 2000_IES, 2010_IES, and 2020_IES) achieved similar results for the four years under study. In all cases, the percentage of pixels classified with an RMSE of less than or equal to 0.025 was between 97 and 98%. In view of the results obtained with respect to the composition of the spectral libraries and the percentage of classified pixels, the results obtained after the classifications made with the optimized monotemporal libraries were analyzed in depth (1990_IES, 2000_IES, 2010_IES, and 2020_IES).

The forest, shrubland, and grassland vegetation classes presented a higher percentage of classified pixels (higher than 20% in all cases). In the case of forest and shrubland classes, more significant differences were observed in terms of the percentage of pixels classified for all years under study. For the grassland class, the percentage was similar for the four years under study: between 26% and 30% of the pixels. The rock and bare soil class did not exceed 20% of classified pixels, whereas the water class did not reach 2% in any of the years under study.

With respect to the type of habitats identified in the study area for each of the analyzed Landsat images (see

Table A1. Table of results obtained from the MESMA fraction images after normalization of the shade showing the number and proportion of pixels for each class model (W = water; F = forest; S = shrubland; G = grassland; R = rock and bare soil).

Class Model	1990		2000		2010		2020
	Pixels	Pixel (%)	Pixels	Pixel (%)	Pixels	Pixel (%)	Pixels	Pixel (%)
W	28,864	0.70	36,038	0.88	31,290	0.76	16,821	0.41
F	93,790	2.29	218,388	5.33	113,279	2.76	521,616	12.72
S	1,012,843	24.71	882,808	21.54	1,268,155	30.94	741,478	18.09
G	523,518	12.77	449,791	10.97	1,189,866	29.03	667,869	16.29
R	124,014	3.03	82,457	2.01	104,392	2.55	49,638	1.21
W,F	33,403	0.81	14,799	0.36	5678	0.14	164	0.00
W,S	27,407	0.67	6944	0.17	4485	0.11	516	0.01
W,G	29,347	0.72	6750	0.16	12,811	0.31	200	0.00
W,R	1945	0.05	1521	0.04	4251	0.10	677	0.02
F,S	494,704	12.07	482,799	11.78	248,272	6.06	631,954	15.42
F,G	261,234	6.37	465,207	11.35	223,935	5.46	325,663	7.94
F,R	235,712	5.75	317,045	7.73	27,639	0.67	45,591	1.11
S,G	392,637	9.58	368,166	8.98	554,261	13.52	585,146	14.27
S,R	456,256	11.13	329,168	8.03	96,066	2.34	164,581	4.01
G,R	332,023	8.10	338,040	8.25	189,882	4.63	253,017	6.17
W,F,S	-	-	-	-	4	0.00	3	0.00
W,F,G	22	0.00	5	0.00	3	0.00	-	-
W,F,R	505	0.01	305	0.01	80	0.00	18	0.00
W,S,G	25	0.00	3	0.00	1	0.00	-	-
W,S,R	20	0.00	16	0.00	6	0.00	2	0.00
W,G,R	39	0.00	22	0.00	33	0.00	4	0.00
F,S,G	7	0.00	15	0.00	1	0.00	26	0.00
F,S,R	31	0.00	32	0.00	2	0.00	4	0.00
F,G,R	171	0.00	164	0.00	36	0.00	272	0.01
S,G,R	259	0.01	340	0.01	63	0.00	308	0.01
TOTAL	4,099,390	100	4,099,390	100	4,099,390	100	4,099,390	100

in Appendix A), the mixture of endmembers of forest and grassland classes modelled the highest percentage of pixels in the four images under study, with 12.35% of the pixels for the 1990 image, 17.78% for the 2000 image, 10.41% for the 2010 image, and 7.94% for the 2020 image.

With respect to the mixture of three types of endmembers, the most frequent was that of forest, grassland, and rock and bare soil classes in the case of the images for the years 1990, 2000, and 2010 (0.02%, 0.04%, and 0.17% of pixels, respectively) and water, forest, and rock and bare soil for the image for the year 2020 (0.01%). The habitats identified as dominant in the study area were shrubland and grassland in the case of the images from the years 1990, 2000, and 2010, and forest and shrubland in the case of the image from 2020.

3.3. Validation Results

The optimized monotemporal libraries provided a lower average value for the RMSE after linear regressions between the reference data and the fraction image data (

Table 4. RMSE average values and standard deviation (SD) obtained of the linear regressions for the different optimized libraries and for each type of habitat and year analyzed.

LIBRARIES		Monotemporal_IES	L1_IES	L2_IES	L3_IES	L4_IES
		0.303 ± 0.089	0.307 ± 0.094	0.318 ± 0.096	0.327 ± 0.107	0.352 ± 0.109
CLASSES		Water	Forest	Shrubland	Grassland	Rock and Bare Soil
YEAR	1990	0.12 ± 0.02	0.30 ± 0.03	0.49 ± 0.04	0.43 ± 0.04	0.26 ± 0.01
	2000	0.26 ± 0.01	0.23 ± 0.03	0.41 ± 0.06	0.43 ± 0.07	0.33 ± 0.02
	2010	0.22 ± 0.02	0.31 ± 0.06	0.45 ± 0.05	0.42 ± 0.05	0.29 ± 0.01
	2020	0.25 ± 0.01	0.24 ± 0.03	0.35 ± 0.06	0.39 ± 0.06	0.26 ± 0.01

). The highest value was obtained for the smallest multitemporal spectral library with equitable distribution (L4_IES).

Furthermore, differences were observed between the mean values of RMSE obtained according to the type of habitat analyzed, with a higher value for the shrubland and grassland classes. The lowest average value of RMSE was obtained for the water class in the classification made with the image of the year 1990. In all other cases, the average values for the RMSE were greater than 0.20. The rock and bare soil and forest had similar mean values between 0.20 and 0.35. The highest mean values of the RMSE were obtained for the shrubland and grassland, exceeding 0.40 in the images of the first three decades.

The mixed linear model applied to analyze the influence of the size of the spectral libraries, as well as their temporality, showed that the size of the library had significant positive effect (p-value = 0.01) on the MESMA classifications. With respect to the temporality factor, the results of the ANOVA of the mixed models showed that its influence was also significant (p-value = 0.01). The RMSE value predicted using the monotemporal libraries was lower than in the case of the predictions made using multitemporal libraries (Figure 7).

On the other hand, no significant influence of the composition of the multitemporal libraries was observed in the obtained results, with differences mainly based on their temporality.

The results presented above show that the monotemporal libraries can be considered the optimal libraries for classification. The monotemporal libraries resulted in the lowest average RMSE values, and the applied statistical models reflected an influence of temporality on the results of the MESMA classifications.

4. Discussion

In this work we analyzed the influence of spectral libraries on the accuracy of a spectral mixture algorithm. Specifically, we analyzed the influence of spectral library size (large vs. small), composition (equitable vs. inequitable distribution of spectral signatures from Landsat TM and Landsat OLI sensors), and temporality (monotemporal vs. multitemporal) on the accuracy of subsequent MESMA predictions. The results of the present work highlight the need to consider the elaboration of the optimal spectral libraries as an essential step to obtain accurate and realistic results.

According to the reported results, larger libraries give rise to more accurate classifications, possibly because their optimal size is determined by the inclusion of a sufficient number of spectral signatures to include the existing spectral variability in a study area [9]. In some studies, libraries have been augmented to account for within-class spectral variability and therefore to improve MESMA analysis [15]. In other cases, pruning techniques have enabled the inclusion of the most representative endmembers of a certain study area in the libraries, reducing their size and heterogeneity and enhancing MESMA performance [37].

The results of the present work show that larger spectral libraries enabled a greater number of MESMA models, resulting in a higher percentage of classified pixels in most cases. However, differences were observed across analyzed years with the same spectral library. The main results obtained in this work showed that a considerable number of spectral signatures were included in the spectral libraries to represent classes such as shrubland and grassland, with even more spectral signatures the optimized multitemporal spectral libraries (L1_IES, L2_IES, and L3_IES). The influence of the high spectral variability within such vegetation classes could affect the selection of a greater number of spectral signatures for each class and could confirm the difficultly of identifying spectral signatures that characterize such vegetation classes [33].

Despite the positive effect of including a variety of endmembers representing within-class spectral variability [38], the use of spectrally similar signatures or the inclusion of many endmembers for the same class has been related to a decrease in the accuracy of classifications [31,39]. In this context, we corroborated the suitability of the IES technique for the optimization of spectral libraries by reducing the size of the starting spectral libraries, obtaining satisfactory results in terms of the Kappa coefficient (greater than 0.9 in all cases) compared to the values of said coefficient obtained in previous studies [16,40].

With respect to the effect of the composition of the spectral libraries, no clear influence of the composition of the libraries was observed in the results. Some previous studies concluded that mixed spectral libraries achieved a similar classification accuracy compared to libraries with data from a single source; however, results with higher precision were also observed when using mixed spectral libraries or spectral libraries with data from a single source [31]. This phenomenon suggests that data from similar sensors can be integrated for such analysis, which would be useful in cases in which data from a single source could be incomplete, such as with Landsat sensors and time series data.

The number of classes of habitats defined in the present work was greater than that found in the analyzed literature, with a total of five classes in the study area; however, models of a maximum of three types of habitats plus shade were applied for the classification of the images, as it has been observed that high similarity between classes can decrease the precision of classifications, resulting in improved when modelling the images with models of two or three classes plus the shade fraction [41]. Therefore, MESMA is limited to modeling fewer classes, even with highly optimized libraries.

The accuracy achieved for classes such as shrubland and grassland (RMSE ≥ 0.35) for the various fraction images was lower than that achieved for the other classes. In addition, the highest accuracy was achieved for the water class was greater, as the spectral signature of this element differs considerably from that of the other elements, with minimal variability, enabling improved identification [42]. The validation of the obtained reflected a significant influence of the size of the spectral libraries in the classifications, as well as their temporality.

5. Conclusions

Spectral unmixing techniques, such as MESMA, allow for the identification of habitat fractions at the subpixel level. However, a necessary step to accurately implement these methods is the generation of optimal spectral libraries.

In the present study, we highlighted the relevance of the size and composition of spectral libraries in determining the performance of spectral mixture analyses. In general, we found that appropriately pruned large spectral libraries led to more accurate MESMA fraction images than small spectral libraries. In addition, we demonstrated the possibility of combining spectral signatures from multiple years and Landsat sensors into a single multitemporal library for MESMA in Landsat time series. However, we demonstrated that monotemporal spectral libraries performed better than multitemporal spectral libraries, suggesting year-by-year analyses to obtain the most accurate results.

Finally, we underline that MESMA methods are subject to limitations when more than three classes of habitats coexist within a pixel. Thus, we encourage further research to address this shortcoming, as the advancement in spectral unmixing techniques is essential to the characterization of actual fine-grained landscapes, as well as to the characterization of past landscapes when pixel sizes of available satellite images are generally equal to or coarser than 30 m.