**Modeling Short-Term Landscape Modification and Sedimentary Budget Induced by Dam Removal: Insights from LEM Application**

#### **Dario Gioia 1,\* and Marcello Schiattarella <sup>2</sup>**


Received: 3 October 2020; Accepted: 27 October 2020; Published: 30 October 2020

**Abstract:** Simulation scenarios of sediment flux variation and topographic changes due to dam removal have been investigated in a reservoir catchment of the axial zone of southern Italy through the application of a landscape evolution model (i.e.,: the Caesar–Lisflood landscape evolution models, LEM). LEM simulation highlights that the abrupt change in base level due to dam removal induces a significant increase in erosion ability of main channels and a strong incision of the reservoir infill. Analysis of the sediment dynamics resulting from the dam removal highlights a significant increase of the total eroded volumes in the post dam scenario of a factor higher than 4. Model results also predict a strong modification of the longitudinal profile of main channels, which promoted fluvial incision upstream and downstream of the former reservoir area. Such a geomorphic response is in agreement with previous analysis of the fluvial system short-term response induced by base-level lowering, thus demonstrating the reliability of LEM-based analysis for solving open problems in applied geomorphology such as perturbations and short-term landscape modification natural processes or human impact.

**Keywords:** landscape evolution; soil erosion; DEM analysis; applied geomorphology; dam removal; base-level lowering; southern Italy

#### **1. Introduction**

Base-level variation has a significant impact on a geomorphological system with severe changes in channel incision rate, sediment flux, and spatial distribution of geomorphological processes [1–3]. A fast transition from endorheic (i.e.: centripetal drainage or closed basin) to exorheic drainage is one of the most relevant cases of disequilibrium of a landscape, because it promotes a non-linear response of the fluvial systems and a complex spatial and temporal response of river incision and sediment flux (see for example [4]). Several works have investigated the long-term response of the drainage network to a transition from endorheic to exorheic conditions due to complex climateor tectonic-driven processes such as sediment overfilling, headward stream erosion, and threshold incision or fluvial capture [5–8]. This kind of analysis is largely based on morphotectonic studies and related qualitative reconstruction of past stages of landscape evolution [9–11]. Dam construction and/or removal is one of the human-induced perturbations of the fluvial net with a stronger impact on the geomorphic system [12–14]. Many works have been focused on the analysis of the effects of this kind of disturbance on the fluvial network and sediment flux, highlighting a typical response that is strongly controlled by the post-dam river longitudinal profile. In fact, upstream knickpoint retreat generally promoted incision of both upland areas and reservoir infill, but several works have

demonstrated that local factors such as rates of knickpoint migration [15,16], bedrock erodibility [17,18], grain-size and texture of reservoir sediment [19] and width of the reservoir [17,19] can drive a complex response of river processes and promote a high spatial and temporal variability of erosion and deposition. Most of these studies are based on field and remote sensing data at limited spatial and temporal scales [17,19,20] or on the application of 1-D models [21,22]. For example, extensive geomorphic analyses of dam removal scenarios have been conducted in the USA through: (i) the direct measurements of the post-dam modification of morpho-sedimentary features supported by one-dimensional hydrodynamic modeling [17]; or (ii) the application of channel evolution models such as CONCEPTS and DREAM [21,22]. Such models have a significant limitation that they are not able to fully capture the spatial (i.e., lateral) distribution of river processes, channel modification, and sediment flux.

In the last years, landscape evolution models (LEMs) have been extensively used to simulate short- and long-term geomorphic scenarios of sediment flux and topographic modification in different natural environments [23–27]. These studies have demonstrated that LEMs represent a powerful tool to predict morpho-sedimentary adjustments related to changes in land use, climate setting, and base levels [24,28–31]. Recently, the Caesar–Lisflood LEM was used to investigate the short-term (i.e., at a decadal scale) pattern and rates of geomorphic changes and associated sediment flux induced by removal of multiple dams in the middle reach of Kaja River, Austria [32]. The authors demonstrated the usefulness of the LEM to predict the 2-D geomorphological and sedimentary effect of the abrupt base-level fall related to dam removal. This kind of investigation overcomes the clear limitation of 1-D or empirical models and can provide a significant opportunity to modify different controlling factors such as climate setting, land use, vegetation growth, and sediment features. The main limitation of the application of such models is that a reliable calibration of a LEM is really complicated and several well-known reasons such as the difficulty of both the selection of model parameters and the absence of extensive validation of the prediction ability of the model results in natural environments could represent a significant limitation in its application in complex natural landscapes [25,33]. Thus, a test area where a LEM is already calibrated can provide a reliable opportunity to simulate the impact of dam removal on the morpho-sedimentary processes and rates. In this paper, we exploited this opportunity and investigated the geomorphic changes induced by dam removal in a small artificial reservoir of southern Italy (Figure 1), where the prediction ability of the Caesar–Lisflood LEM has been already tested [34]. More specifically, the calibration of the model has been recently carried out through the comparison between the simulation results and a direct estimation of sedimentation rates in the reservoir, demonstrating a good prediction ability of the model [34]. Application of the dam removal scenario in the catchment allowed us to investigate the short-term (i.e., at a decadal scale) response of the geomorphic system to the abrupt change in base level as well as the role of several local factors (i.e., infill thickness and geometry, the sediment and bedrock features) on the changes in sediment flux and channel pattern.

#### **2. Study Area**

The study area is located in an upland area of the Ofanto basin, a Pliocene-Quaternary tectonic depression of the southern Apennines, Italy (Figure 1). It includes the catchment of the Ficocchia Torrente stream, a dextral tributary of low hierarchical order of the Ofanto River. The Ficocchia stream is dammed in its lower reach by an earth dam (i.e., the Saetta dam), which was constructed by EIPLI (Agency for the Development of the Irrigation and Agricultural Transformation, Ministry of Agricultural, Food and Forestry Policies) between 1988 and 1989.

The Ofanto River cuts a large E-W-trending intermontane tectonic depression, Pliocene to Quaternary in age. The infill of the basin is intensively deformed by Pliocene-Pleistocene folding and faulting related to the younger stages of evolution of the southern Apennines chain [35] and covers a wide (i.e., about 350 kmq) and elongated area along the Ofanto River valley [36,37]. Pliocene-Pleistocene deposits of the Ofanto basin are composed of clay, sandstones and conglomerates, which unconformably overlay poly-deformed units of limestone, shale and sandstone that belong to the Irpinian and Lagonegro tectonic units of Cretaceous to Miocene ages [36].

**Figure 1.** Geological map of the study area (modified from [36]). Legend: (**1**) Clay of lacustrine environment (lac, Holocene) (**2**) Landslide deposits (lan, Holocene) (**3**) Holistolits made by decametric blocks of limestone (pa, Upper Miocene); (**4**) Silt and marly clay (CVT2, Upper Miocene); (**5**) Coarse- to medium-grained sandstone with rare intercalation of lens of polygenic conglomerate (CVT1, Upper Miocene); (**6**) Calcareous breccia and grey shale (FYRa, Lower Cretaceous-Oligocene); (**7**) Alternance of chert, marly clay, calcarenites and calcareous breccia (FYR1, Lower Cretaceous-Oligocene); (**8**) Light-grey and greenish shale with intercalation of marls and limestone (FYG, Lower Cretaceous); (**9**) Alternance of calcarenite, calcilutite and varicoloured clay (FMS, Upper Cretaceous-Eocene); (**10**) Varicoloured clay (AVF, Lower Cretaceous); (**11**) High-angle fault (dashed if uncertain); (**12**) Thrust (dashed if uncertain); (**13**) Stratigraphic contact. (**A**) Geographical location of the study area.

The reservoir catchment is carved in a Cretaceous to Miocene bedrock made by tectonic domains belonging to Lagonegro units, Sicilide calcareous-clay succession and flysch deposits of Miocene syntectonic basins (Figure 1, see also [36]). The study area also features widespread outcroppings of upper Miocene deposits belonging to Castelvetere Formation. These deposits are composed of light brown sandstone with intercalation of conglomerate (CVT1, Figure 1) passing upward to silt and marly clay (CVT2, Figure 1) containing large blocks of olistoliths (pa, Figure 1). These deposits mainly crop out in the south-western sector of the study area and unconformably overlies the Lagonegro tectonic units, which form an elongated belt in the north-eastern sectors of the study area. The outcropping deposits of the Lagonegro units are mainly constituted by lower–middle Cretaceous marls and grey shales with calcarenites and calcirudites and calcareous breccia (Flysch Rosso Formation, FYRa and

FYR1 in Figure 1). Lower Cretaceous varicoloured clay (AVF) with thin intercalations of calcirudites and calcarenites (FMS) crops out in the northern sector of the studied area (Figure 1).

Heterogeneous landslide deposits and fine-grained reservoir sediments are the youngest lithological units of the study area.

The geomorphological pattern of the study area is strongly influenced by the Pliocene-Quaternary tectonic evolution and relief growth of this sector of the chain, which have also controlled the lithological features and spatial distribution of deposits. The landscape is located at an altitude ranging from 942 to 1242 m a.s.l. and is dominated by E–W trending ridges and thrust sheets, which are mainly carved in Cretaceous-to-Miocene pelagic deposits. These structural landforms are deeply cut by the drainage network of the Ficocchia stream (Figure 2).

**Figure 2.** (**A**) DEM of the pre-dam removal landscape and drainage network of the study area. Hierarchization follows the Strahler's scheme. Numbering of the catchments is shown in the frame to the left. (**B**) Land-use map. Legend: (**1**) Anthropic surfaces and roads; (**2**) Arable lands; (**3**) Sclerophyllous vegetation; (**4**) Broad-leaved and mixed forests; (**5**) Natural grasslands; (**6**) Water courses and water bodies. (**C**) Isopachs of soil thickness derived by a GIS-based interpolation of the results of a field-survey analysis. Modified after [34].

The drainage basin of the artificial reservoir is formed by three small catchments of low hierarchical order showing a well-organized drainage network with a sub-dendritic pattern. The watershed of the southern sector of the reservoir catchment is featured by a sub-circular shape and runs on low-relief erosional land-surfaces [34,37]. In this sector, headwater channels exhibit a higher gradient than the channels located in lower reaches (Figure 2). The main geomorphological processes of the study area are channel incision and fluvial erosion, which also represent the main influence factor of the sediment yield and flux. Slope processes and landsliding phenomena related to minor shallow mass movement processes and small earth-flows are located in the eastern sector of the studied catchment, where clay-rich deposits crop out.

Regional climate data of the study area was derived from the statistical analysis of a weather station located about 20 km to the west of the study area (daily rainfall record of the Pescopagano weather station, period: 1951–2015, [38]). Mean annual rainfall of the study area in the last half-century is slightly higher than 1000 mm year−1. The climate setting is featured by dry summers and cold winters with a maximum of rainfall in autumn. In the last decades, a general trend of decrease in total annual rainfall was observed. This trend is associated with a decrease in rainy days and an increase of multi-day extreme rainfall events, mainly in autumn and spring. Shorter-term rainfall record (period: 1994–2016) of the hydro-meteorological station of the study area highlights a similar climate trend: mean annual precipitation is about 900 mm and rainfall peaks are mainly observed in the autumn and winter seasons, with rainfall maxima of 150–160 mm [34].

According to the classification of the III level of the Corine Land Cover project [39], a land-use map of the catchment was prepared from a revision of literature data and original investigation based on photointerpretation of aerial and satellite images. The landscape is dominated by semi-natural areas with natural grasslands and sclerophyllous vegetation that cover about 85% of the total area. Other sectors of the study area are classified as agricultural and urban areas.

#### **3. Materials and Methods**

#### *3.1. Caesar–Lisflood LEM: Model Description and Calibration*

Caesar–Lisflood is a second-generation LEM that can be used to model short- and long-term topographic changes and sediment flux of complex natural landscapes. In the last decade, the model was extensively used for the analysis at different temporal and spatial scales of issues of applied geomorphology and hydrology. As a matter of fact, there are at least 60 published papers dealing with model applications [40]. A rigorous and complete description is beyond the scope of this paper and is reported in [41].

The Caesar–Lisflood LEM uses a hydrological model to generate spatially distributed runoff, which is generated on a DEM to estimate flow depths and velocities [41]. Such data are used to model topographic changes and to assess topographic modification (i.e., erosion and deposition) within an active layer with a pre-defined grain-size of deposits. The modified DEM becomes the starting point for the next time step of the simulation. The user can select a catchment mode with no internal influxes other than rainfall or a reach mode, where the flow enters through the main river. Input data of the model are: DEM, grain-size features, rainfall, depth of bedrock and Manning coefficient values [24,41].

Rainfall precipitation represents the input to derivate catchment runoff, which in turn drives topographic changes (i.e., erosion and deposition) and controls fluvial and slope processes for the modeled time step [28]. Depths and velocities of the flow are estimated from flow discharges between raster cells using the Manning's equation and are then used to model both the sediment transport and erosion/deposition processes. Caesar–Lisflood estimates sediment transport in relation to nine grain-size fractions, which are selected by use. Sediment can be transported as bed load or as suspended load. Soil creep and landslide processes can be also included in the simulation [41,42] using a critical slope angle threshold. This condition allows the re-distribution of landslide and soil creep deposits from slopes to the fluvial system. Elevations and grain size features of the cells are updated according to the estimation of erosion/deposition and slope process. Model outputs are: (i) elevation changes; (ii) flow discharges; (iii) sediment fluxes at the outlet over time [24,28,41–43]. The model is able to reconstruct topographic changes at a sub-metric scale. Input data and model parameters represent the key to derive a robust predictive model, although a recent work by [44] demonstrated that results of

the Caesar–Lisflood LEM are influenced only by a limited number of parameters. More specifically, the authors suggest that spatial distribution of Manning coefficient and the selection of the sediment transport formula are the main sensitive parameters of the model results. For this reason, we have carefully selected these parameters: for example, Einstein's sediment transport formula [45] was reconstructed using laboratory tests on (predominantly) sand-based channels and we have here chosen this equation according to the prevalent grain-size deposits of the study area. A 5-m DEM was used for the simulation. Moreover, the other input parameters were accurately derived from lithological, climate, and land-use features of the study area. To introduce a representative framework of the study area, we have preliminarily reconstructed the bedrock depth through a field-based analysis of the lithological units. To this aim, a detailed lithological map has been drawn by review of previous works and new geological surveys. Moreover, field surveys of the soil depth have been performed for the different lithological units of the study area (Figure 2, [34]). This approach allowed us to infer the soil thickness and bedrock depth for each lithological. This kind of information was summarized in a map showing the isopachs of the mean soil thickness (Figure 2c). The map has been converted in a raster grid and was introduced in the model in order to set the depth of bedrock and grain-size distribution within the active layer where erosion and deposition processes are estimated by LEM. The layer of the bedrock has been derived in a GIS environment using a map algebra tool based on a subtraction between the soil thickness map and the DEM.

Hourly rainfall dataset was derived from a weather station located near the catchment outlet. Rainfall record covers a time interval of about 22 years, and it was used also for the damremoval scenario.

Input data and model parameters are summarized in Table 1 (see also [34]). It is worthy of note that the estimation of the sedimentary budget base on the LEM simulation was already tested through a comparison between the total amount of eroded sediment volumes coming from the model results and direct measurement of the short-term (i.e., about 20 years) sediment storage within the artificial reservoir. The good accordance between the model results and direct measurements [34] suggests that the model can be used to simulate the dam removal scenario in a consistent way.



Manning coefficient values were assigned (Table 2) according to a revision of the values proposed in previous works (see for example [46,47]) whereas the spatial distribution of the Manning coefficient follows the boundaries of the land use map of the study area (Figure 2b).


**Table 2.** Values of the Manning coefficient.

#### *3.2. Dam-Removal Scenario*

Two different scenarios were simulated in the study area (Table 3): the first one is the pre-dam removal scenario, where the modification of the initial topography (PreDR-T0, Table 3) has been simulated over a time interval of 20 years (Pre-DR-T20, Table 3). Such a scenario represents the present-day landscape and the simulation provided a reconstruction of the sediment flux in the reservoir and the short-term topographic changes of the fluvial and slope systems. Input parameters and boundary conditions are summarized in Table 1.

**Table 3.** Description of the simulation scenarios: the first scenario refers to the present-day geomorphological setting whereas scenario 2 infers the geomorphic changes induced by the dam removal.


The post-dam removal scenario has been developed removing the dam body to the initial DEM. More specifically, the DEM was created by subtracting to the original topography the dam height (PostDR-T0, Table 3). Post-dam removal simulation was run for a period of 20 years (PostDR-T20) using the same hourly rainfall data and input parameters of the Pre-DR scenario.

Output DEMs and raster of the topographic changes were analysed in a GIS environment to assess the pattern and rates of sediment flux and channel profile adjustments. Spatial distribution of the erosion/deposition processes, multi-temporal analysis of river longitudinal profiles, and valley topographic changes represent the key data to infer the geomorphic response of the study area to the simulated fall of the base-level.

#### **4. Results**

#### *4.1. Scenario 1, Pre-Dam Removal*

Outputs of the Caesar–Lisflood LEM are analysed in a GIS environment in order to investigate the geomorphic changes of the fluvial system and the spatial and temporal distribution of sediment erosion and deposition and their relationships with lithology, land use, and geomorphological processes. Figure 3 shows the hillshades representing the initial DEMs used for the modeling of the two simulation scenarios. Dam removal promoted an increase of the catchment area of about 0.6 kmq. The post-dam removal catchment includes the steeper reach of the Ficocchia stream, which flows in a deep V-shape valley.

As already described in the previous section, a recent work provided a comparison between LEM-based estimation of sediment flux of the Ficocchia catchment and direct measurements of reservoir sedimentation volumes [34]. This estimation covers a period (i.e., 18 years) similar to the simulation period of this work and highlights good accordance between model results and field-based data. More specifically, LEM-based erosion volumes tend to underestimate the reservoir sedimentation of about 20%.

**Figure 3.** Hillshades representing the initial DEMs used for the modeling of the different simulation scenarios: (**A**) pre-dam removal, initial topography (PreDR-T0); (**B**) post-dam removal, simulation period: 1 year (PostDR-T0).

Simulation results of the pre-dam removal scenario are summarized in a DEM of difference (Figure 4), which shows the topographic changes at a sub-metric scale and sediment flux in the study area after 20 years. Visual inspection of the map allowed us to infer the short-term morpho-sedimentary evolution of mountain catchments and the influence of fluvial and hillslope processes on the geomorphological evolution and sediment delivery. The map highlights that erosion processes occurred along the main channels of the fluvial net whereas deposition is located in the artificial reservoir as a result of the flattening of longitudinal profiles in the lower reaches of the channels (Figure 4).

**Figure 4.** (**A**) Elevation difference map for the pre-dam removal scenario (simulation period: 20 years); (**B**) Numbering of the three sub-basins of the study area.

Model results seem to suggest a minor impact of slope processes on the sediment yield of the studied catchment, although landslide processes and high levels of the reservoir water table have been invoked as possible factors of slight differences between LEM-based erosion volumes and direct estimation of reservoir sedimentation [34]. In any case, modelling results fit well with the spatial distribution of landforms and deposits deriving by field-based geomorphological analysis. In fact, several geomorphological evidences such as V-shape valleys of the main channels and the absence of slope and alluvial deposits along the thalwegs is in accordance with the LEM results and suggest that channel incision upstream of the artificial reservoir is the main geomorphological processes controlling the morpho-sedimentary evolution of the Ficocchia catchment. As a matter of fact, sediment delivery ratio of the artificial reservoir is slightly higher than 0.9, testifying a high level of sediment connectivity of the study area.

Histogram of Figure 5 shows the results of the elevation changes caused by erosion and deposition after a time interval of 20 years. The total amount of erosion deriving by Caesar–Lisflood LEM simulation for the entire period is 84,070 m3, which corresponds to a mean annual erosion volume of 4203 m3/year (Figure 5). Stable areas represent 97.7% of the total area whereas the most representative erosion class has a value ranging from 0.1 to 0.5 m (Figure 5). Analysis of the contribution of each sub-basin to the sediment yield highlights that about 83% of the total amount of erosion volumes comes from catchment n.2 and n.4 (Figure 5).

**Figure 5.** Erosion/deposition classes in the catchment deriving from the analysis of the altitude difference map of Figure 4. In the frame: distribution of the eroded volumes from the three sub-basins of the study area (numbering is shown in Figure 4).

#### *4.2. Scenario 2–Dam Removal*

Post-dam removal scenario (Post-DR, Table 3) was carried out over a 20-year period to perform a direct comparison with the sediment budget related to Scenario 1. Decadal-scale modelling of the geomorphic response to dam-removal predicts a strong modification of the fluvial system, with a general increase of topographic changes than the pre-dam removal scenario (Cfr. Figures 4 and 6). Altitude difference map (Figure 6) clearly highlights a significant increase of erosion processes in both the headwaters and mid sectors of the drainage net, which are mainly related to the higher incision ability of the main channels. In fact, a comparison between the results of the two simulations shows that lateral migration of channels appears to be limited (cf. Figures 4 and 6). LEM predicts the development of a wide floodplain in the flat area of the reservoir with a well-defined incision of the infill by the three main channels of the study area. In addition, analysis of the spatial pattern of erosion and deposition shows the occurrence of small sedimentation areas upstream to the reservoir infill (Figure 6), which contributed to a slight decrease of the sediment connectivity of the entire catchment.

**Figure 6.** (**A**) Landscape evolution models (LEM)-based elevation difference map for the post-dam removal scenario (simulation period: 20 years). (**B**) Numbering of the post-dam removal catchments.

The frequency distribution of the altitude difference map (Figure 7) illustrates a higher erosional ability of the geomorphic system after the lowering of the base level.

**Figure 7.** Statistical distribution of the altitude difference map of Figure 6 (Post-dam removal scenario). In the frame: distribution of the eroded volumes from the three sub-basins of the study area (numbering is shown in Figure 6).

In fact, the model predicts an increase of the total eroded volumes by a factor of 4.4 for the post-dam removal scenario whereas the mean annual erosion volumes increase from 4203 m3/year to 18465 (Figure 7). Caesar–Lisflood simulation infers topographic modifications for a percentage of the total area of about 7.1% (corresponding to an area of 0.74 kmq). In this scenario, erosion classes with a value ranging from 0.1 to 0.5 m represent the maximum of the altitude difference map but the comparison with Scenario 1 shows a relative increase of the erosion classes with a higher value. Analysis of the spatial pattern of erosion for the different sectors of the study area shows that sub-basin 4 and 3 contribute to 45.7% and 18.7% of the total amount of eroded volumes yielded in the upland sectors.

Then, approximately 81.3% of the total amount of sediment eroded comes from sub-basin n.3 and n.4 (i.e., the easternmost and the central one), highlighting a higher erosion ability of these catchments for the dam removal scenario.

#### *4.3. River Profile Analysis and Channel*/*Valley Modifications*

Multitemporal analysis of longitudinal river profiles has been performed using the elevation data derived from the output DEMs for each scenario. The perturbation induced by the dam removal promoted a significant modification of the river profiles (Figure 8), with a complex spatial and temporal pattern of geomorphic adjustments. Moreover, our simulation also highlights the occurrence of 2-D channel modifications, which are mainly related to a lateral shift of the main channels in the higher altitude sectors surrounding the reservoir infill (Figure 8).

After the dam removal, the longitudinal profile evolution is mainly featured by a pronounced channel incision, which mainly occurred upward the flat area of the reservoir (Figure 8). In fact, the fast response of main channels to the base-level lowering is a lengthening of the longitudinal profiles, which promoted incision upstream of the former reservoir area (Figure 8, see the yellow curves and the black ones). Channel 2 adjusts its longitudinal profiles forming a major convex knickzone downstream of the removed dam, which retreats of about 180–200 m during the 20-year simulation period. A minor knickpoint can be observed in the three channels at an altitude of about 955–960 m, showing a lower rate of retreat (i.e., about 1 m/yr).

Cross profiles of Figure 9 furnished additional information about the landscape modification resulting from our simulation. Spatial and temporal evolution of valley/channel modification is complex, with alternating stages of incision and aggradation. Profiles located in the higher altitude sectors of the catchment are featured by meter-scale incision (see for example profile a-a', d-d' and e-e'); in this sector, the highest amount of deepening (i.e., about 3 m) can be observed at profile a-a'. Topographic profiles crossing the reservoir deposits shows a more complex response to dam removal due to a general tendency of channels to cut the infill in the southernmost sectors (Figure 9, see profile b-b' and f-f') and the occurrence of aggradation upward the removed dam (see for example the profile g-g', Figure 9). Pronounced incision phenomena occurred downstream of the dam, with a progressive deepening of the main channels of about 3 m over the 20-year simulation period (Figure 9, profile i-i' and l-l').

**Figure 8.** Map showing the planimetric changes of the main channels from the present-day landscape (initial DEM, Pre-DR-T0 in Table 3) to the final stage of the post-dam removal scenario. To the bottom: comparison of longitudinal river profiles of the three main channels (channel 1, 2 and 3 in the map) deriving from the simulation scenarios. River profile analysis highlights the amount of incision related to the base-level fall as well as the development of pronounced knickpoints in the upper and lower reaches of the main channels. Higher rates of knickpoint retreat are associated with the lower reach of the channel 2.

**Figure 9.** Topographic profiles for the different simulation scenarios (location of the profile is reported in the map) showing the landscape modification resulting from the simulation. Higher rates of fluvial incision occurred in the upper and lower reaches of the catchment (about 3 m over the 20-year simulation period, see for example profile a-a' and i-i'). A deep incision of the reservoir top can be also observed.

#### **5. Discussion**

LEM-based simulation of the geomorphic response of an upland catchment of southern Italian Apennine to dam removal suggests a complex spatial and temporal modification of channels and sediment flux. Model calibration and selection of the appropriate input parameters and boundary conditions are a critical issue for the LEM application in complex landscapes. (Skinner et al., 2018) [44] provided a detailed sensitivity analysis of the Caesar–Lisflood LEM, demonstrating that only some parameters such as sediment transport formula, Manning coefficient, and sediment grain sizes have a significant influence on model results.

Our approach based on a calibration of the model using direct measurements of sedimentation volumes in the artificial reservoir overcomes the issue of the complex selection of an excessive number of input and model parameters of the Caesar–Lisflood LEM, which are frequently ascribed as the main factor of its rare application [33,34,44].

In the study area, these parameters have been here derived by detailed field measurements and calibrated by the evaluation of an independent dataset. The preliminary validation of the short-term prediction ability of the Caesar–Lisflood LEM deriving from a direct comparison between model results and direct measurements of reservoir sedimentation volumes is a strong point of the proposed approach, which suggests the robustness of the model settings and parameterization. In fact, the limited difference between source-sink data confirms the robustness of the model and allowed us to consider our simulation as a robust estimation of topographic changes, geomorphological processes, and sediment flux at a short-term scale (i.e., 20 years) induced by base-level lowering. Such an approach can be effectively integrated with short-term analysis of geomorphic processes based on multitemporal comparison of high-resolution DEMs [48,49] and/or quantitative geomorphological analyses [8,50–52].

Our investigation demonstrated that Caesar–Lisflood LEM has a high potential to explore scenarios of morpho-sedimentary changes in response to different perturbing factors such as base-level fall or climate/land-use changes. The main advantage of the use of the LEM is its ability to reconstruct 2-D topographic modification and sediment flux, which can be investigated through advanced GIS-tools of map algebra and spatial statistics.

The geomorphic response to dam-removal scenario of the Ficocchia catchment mainly consists of a higher erosion ability of the channels in the upper reaches of the fluvial systems and a deep incision in the lowermost sectors of the post-dam removal catchment. As a matter of fact, quantitative analysis of the sediment dynamics resulting from the dam-removal scenario highlights a significant increase of the total eroded volumes in the post dam scenario of a factor higher than four. This significant increase of sediment yield is a common response of the fluvial system to base-level lowering, although our predicted erosion rates of the study area after the dam removal are significantly higher than the estimation coming from previous works (see for example [32]). For example, a recent analysis of the impact of multiple dam removal using Caesar–Lisflood LEM [32] reported a lower amount of sediment delivery induced by such a dam removal scenario.

This kind of response of the fluvial system to base-level lowering is in accordance with reconstruction made by other studies since it is mainly represented by upstream widespread incision associated to knickpoint retreat [17,19,21,32].

Predicted rates of knickpoint retreat reached values lower than those predicted by other researchers [32,53]. For example, a recent LEM-based analysis of multiple removals of narrow dams in Austria predicted annual knickpoint retreat rates ranging from 150 to 300 m [32], which is an order of magnitude higher than our results. As a matter of fact, our analysis highlights a peak of knickpoint retreat for the main channels of about 10 m/yr, see Figure 9).

The presence of a wide flat landscape coinciding with the former reservoir could represent the external factor that drives the complex and peculiar response of the study area. This sector can promote the dysconnectivity between the two steeper segments of the longitudinal profiles located upward and downward the removed dam, which could inhibit a faster retreat of the major knickpoint related to dam removal and promoted the formation of a wide floodplain in the mid sectors of the catchment. This observation confirms the relevant role of local parameters and physiographic setting (i.e., geometry and sedimentary features of the reservoir infill, bedrock and sediment erodibility, channel profile slope, etc.) in controlling the spatial distribution of erosion/deposition processes and the response of the fluvial system to severe perturbations induced by a fast base-level lowering, thus emphasizing the crucial role of landscape evolution models in the reconstruction of complex spatial and temporal changes of erosion and deposition induced by human disturbances.

#### **6. Concluding Remarks**

Simulation scenarios of sediment flux variation and topographic changes due to dam removal have been investigated in a reservoir catchment of southern Italy through the application of a landscape evolution model (i.e.: the Caesar–Lisflood LEM). LEM simulation highlights that the abrupt change in base-level due to dam removal induces a significant increase in erosion ability of main channels and a strong incision of the reservoir infill. Analysis of the sediment dynamics resulting from the dam-removal highlights a significant increase of the total eroded volumes in the post dam scenario of a factor higher than four. Model results also predict a strong modification of the longitudinal profile of main channels induced by dam removal, which promoted deep fluvial incision in their upper and lower reaches.

Results are in agreement with previous analysis of the short-term response of the fluvial system, thus demonstrating the reliability of LEM-based analysis for solving open problems of short-term topographic changes and sedimentary budget induced by natural or human perturbations. The general good accordance between the model results and independent analysis based on field data ([34], see also paragraph 4.1) demonstrated both the general usefulness of the approach for the investigation of human-induced geomorphic disturbance of a landscape and the usefulness of the proposed approach than the application of 1-D simplified models. The main advantage of the use of the Caesar–Lisflood LEM in reconstructing the 2-D spatial and temporal pattern of topographic changes and related geomorphological processes is that the simulation scenarios can be easily compared with field data and historical maps, which can be also useful to explore the source of uncertainties, simplifications, and assumptions of the model.

**Author Contributions:** Conceptualization, D.G. and M.S.; methodology, D.G. and M.S.; software, D.G.; validation, D.G. and M.S.; formal analysis, D.G.; investigation, D.G. and M.S.; resources, D.G. and M.S.; data curation, D.G.; writing—original draft preparation, D.G.: writing—review and editing, D.G. and M.S.; visualization, D.G. and M.S.; supervision, D.G. and M.S.; project administration, D.G. and M.S.; funding acquisition, D.G. and M.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** We would like to thank the three anonymous referees for their helpful suggestions and constructive reviews. The first author would also like to thank Maurizio Lazzari (ISPC-CNR) for support in fieldwork and interesting discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Comparison of Satellite and Drone-Based Images at Two Spatial Scales to Evaluate Vegetation Regeneration after Post-Fire Treatments in a Mediterranean Forest**

**Jose Luis Martinez 1, Manuel Esteban Lucas-Borja 1, Pedro Antonio Plaza-Alvarez 1, Pietro Denisi 2, Miguel Angel Moreno 3, David Hernández 3, Javier González-Romero <sup>1</sup> and Demetrio Antonio Zema 2,\***


**Citation:** Martinez, J.L.; Lucas-Borja, M.E.; Plaza-Alvarez, P.A.; Denisi, P.; Moreno, M.A.; Hernández, D.; González-Romero, J.; Zema, D.A. Comparison of Satellite and Drone-Based Images at Two Spatial Scales to Evaluate Vegetation Regeneration after Post-Fire Treatments in a Mediterranean Forest. *Appl. Sci.* **2021**, *11*, 5423. https:// doi.org/10.3390/app11125423

Academic Editors: Dario Gioia and Maria Danese

Received: 16 March 2021 Accepted: 8 June 2021 Published: 10 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** The evaluation of vegetation cover after post-fire treatments of burned lands is important for forest managers to restore soil quality and plant biodiversity in burned ecosystems. Unfortunately, this evaluation may be time consuming and expensive, requiring much fieldwork for surveys. The use of remote sensing, which makes these evaluation activities quicker and easier, have rarely been carried out in the Mediterranean forests, subjected to wildfire and post-fire stabilization techniques. To fill this gap, this study evaluates the feasibility of satellite (using LANDSAT8 images) and drone surveys to evaluate changes in vegetation cover and composition after wildfire and two hillslope stabilization treatments (log erosion barriers, LEBs, and contour-felled log debris, CFDs) in a forest of Central Eastern Spain. Surveys by drone were able to detect the variability of vegetation cover among burned and unburned areas through the Visible Atmospherically Resistant Index (VARI), but gave unrealistic results when the effectiveness of a post-fire treatment must be evaluated. LANDSAT8 images may be instead misleading to evaluate the changes in land cover after wildfire and post-fire treatments, due to the lack of correlation between VARI and vegetation cover. The spatial analysis has shown that: (i) the post-fire restoration strategy of landscape managers that have prioritized steeper slopes for treatments was successful; (ii) vegetation growth, at least in the experimental conditions, played a limited influence on soil surface conditions, since no significant increases in terrain roughness were detected in treated areas.

**Keywords:** VARI; log erosion barrier; contour-felled log debris; land restoration; vegetation cover

### **1. Introduction**

Wildfires can negatively affect soil fertility, biodiversity, land resources, global warming, and human assets; however, positive environmental effects are also recognized, such as increased forest regeneration and nutrient recycling [1]. The Mediterranean ecosystems are adapted to fire disturbance, but the increasing recurrence and high severity of wildfires generate high soil loss and reduce the ability of vegetation to recover [2].

In the Mediterranean Basin, wildfires burn large forest areas not only during summer, but also in mid seasons, when heavy rainstorms may occur [3,4]. For instance, in Spain, even though the number of wildfires has followed a decreasing trend in the last decades [5], wildland is still severely affected by forest wildfires in summer. In the last 10 years, more than 3000 km<sup>2</sup> of forests have burned. Only in 2018, over 7000 wildfires burned about 300 km2 of forests [5].

The hydrological and ecological effects of large wildfires (floods, landslides and mudslides, erosion, loss of biodiversity, destruction of fauna, etc.) can be devastating [6]. Among these effects, soil erosion is presumably the most severe consequence of forest fires, since it threatens water resources, infrastructures, and populations inside and out of burned areas [7]. For burned soils, erosion after wildfire is much higher compared to unburned, since fire changes several physical and chemical properties of soil [1,7]. For instance, erosion in burned forest areas increases by even 30 times compared to the natural values recorder in unburned zones [8]. To reduce, under tolerable limits, soil erosion, urgent restoration actions are implemented immediately after the wildfire to prevent secondary post-fire damage, whereas long-term restoration can be used for regenerating plants and promoting the evolution of forest structure [9].

Management actions in fire-affected forests have typically focused on soil structure and plant recovery. The increase in vegetal cover through plant regeneration has been identified by some studies in Mediterranean regions as the main driver in post-fire management of forests [10]. Among the post-fire restoration actions, log erosion barriers (hereinafter indicated as "LEB") or contour-felled log debris (CFD) are widespread measures to control the erosive processes (mainly due to runoff and sediment flows) in forest ecosystems [11]. More specifically, LEBs are built by felling burned trees, while CFDs consist of branch and small-felling burned trees, and both are laid on the ground along the slope contour. Many studies have assessed the effectiveness of hillslope stabilization treatments on soil hydrological response (namely runoff and erosion) (e.g., [12–15] as well as reductions in CO2 emissions and carbon sequestration [16]. In some cases, the impact of these actions in reducing runoff and trapping sediments is limited to the less intense rainfall events (e.g., [11]). In other cases, their effects are scarce or negligible, due to inadequate constructions or deficient design [13]. Some studies have even shown a negative impact, in terms of a higher percentage of bare and stony soil, especially in sunny areas subjected to different forest treatments, with subsequent increase in soil erosion [17]. It is evident that the effects of post-fire hillslope treatments are still uncertain, and the most suitable technique has not been completely identified [18].

A monitoring activity of post-fire management actions in burned forests in the Mediterranean environment may give forest managers and landscape planners insight about their effectiveness in reducing soil loss and accelerating the vegetation recovery. This activity is essential, since public administrations make large efforts and spend many resources in restoration actions that usually include economic, cultural and even landscape components [19,20]. However, the techniques of field monitoring are generally expensive and time-consuming [21]. The remote sensing techniques (using satellites and drones) are new tools to make land surveys quicker and cheaper [22]. Remote sensing is usually applied in fire-affected areas for fire risk and fuel mapping, active fire detection, burned area estimates, burn severity assessment, and post-fire monitoring vegetation recovery [23]. The remote sensing techniques using satellite images are able to evaluate post-fire regeneration with increasingly improved spatial, temporal and radiometric resolutions [24]. In the last decade, the appearance of unmanned aerial vehicles has increased the use of remote sensing application to agro-forestry [25,26]. UAV technology is promising for monitoring vegetation regrowth, since the spatial resolution and temporal intervals of surveys are not dependent on satellite orbits [23]. Despite their versatility and low cost, UAVs are not widely used to survey wide areas due to legal and technical limitations (e.g., autonomy, payload capacity) [27]. The assessment of wildfire effects and post-fire regeneration using remote sensing has been done through the calculation of spectral vegetation indices [28]. These indices are calculated using the electromagnetic wave reflectance data of vegetation using passive sensors [29]. Several studies have successfully developed and applied vegetation indices (e.g., the Normalized Difference Vegetation Index, NDVI; the Composite Burn Index, CBI; the Differenced Normalized Burn Ratio, DNBR, and the Visible Atmospherically Resistant Index, VARI [24]. In this regard, Gitas et al. [25] published a comprehensive review about remote sensing of post-fire vegetation recovery and highlighted the important role of remote sensing in the future as well as the need to perform more studies in Mediterranean areas. Bhagat et al. [26] reported that UAVs provide data at fine resolution with the desired temporal resolution, which make them cost-effective and efficient data collectors. Regarding satellite images, Sirin et al. [29] was able to compare different multispectral satellite data to assess vegetation cover in abandoned lands and rewetted peatlands. Bright et al. [30] analyzed post-wildfire recovery in different coniferous forest types in North America using the Normalized Burn Ratio (NBR) derived from LandTrendr images. The same index was used in the same environment by Lentile et al. [31] to monitor post-fire burn severity and vegetation response following eight large wildfires. In the Mediterranean Basin, post-fire vegetation recovery was assessed by Bastos et al. [32] in Portugal (using the VEGETATION sensor), Lasaponara et al. [33] in Italy (using satellite MODIS data), and Polycronaki et al. [34] in Greece (using optical and SAR data [35].

Regarding VARI, several applications have been carried out in different research fields. For example, Schneider et al. [36] showed that VARI outperforms other indexes in distinguishing historical wildfire data in southern California. In the same environment, Stow et al. [37] demonstrated better performances of VARI compared to another normalized index (Normalized Difference Water Index, NDWI) in monitoring chaparral moisture content. Munoz et al. [38] found that VARI was more effective in estimating the fraction of vegetation cover and recognizing the different land use compared to NDVI index, observing a standard error lower than 8%. From these studies, the feasibility of VARI to estimate variable levels of land cover emerges. However, applications of this index for evaluation of vegetation regeneration in fire-affected areas are very scarce. Only Larrinaga and Brotons [39] calculated four greenness indices (among them VARI) to analyze post-fire regeneration of Mediterranean forests. Christakopoulos et al. [40] proposed a comparative evaluation of restoration practices using remote sensing and GIS in naturally-regenerating and reforested areas of Greece. No research is available about the estimation of vegetation recovery after hillslope stabilization techniques using VARI applied to satellite or UAV surveys. This leaves the usability of these methods in post-fire management of forests not well understood, to date.

To fill this gap, the current study evaluates the ability of VARI, estimated by both satellite and UAV images, to quantify the vegetation recovery after post-fire treatments (LEBs and CFDs) in a Mediterranean forest of Central Eastern Spain. Possible correlations between the vegetation cover measured in field plots and VARI values estimated from satellite (LANDSAT8) or UAV images, are identified. Processing of this information at two spatial scales (catchment and hillslope) allows validation of one or both methods.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The investigation was carried out in the Sierra de los Donceles (municipality of Hellín, Castilla-La Mancha region, southeast Spain) (Figure 1a,b). This mountain forest lies in the most northeastern part of the Penibetic system with a southeast-northwest aspect. The landscape is patchy with steep slopes (sometimes over 25–30%), predominantly exposed to south-southeast, and ranges between 304 and 808 m above sea level with an average altitude of 506 m.

In this forest, one catchment (hereinafter "catchment area" or "catchment scale") and one hillslope ("hillslope area" or "hillslope scale") were selected, both subject to burning and post-fire restoration actions (Figure 1c,d). In these areas, plots for survey of vegetation were installed, as detailed in Section 2.2, and remote sensing-based analysis (using satellite or drone for the catchment and hillslope scale, respectively).

**Figure 1.** Location (**a**) and aerial map (**b**) of the study area, catchment area (**c**) and hillslope area (**d**) in Sierra de los Donceles forest (Castilla La Mancha, Spain).

The climate of the region is semi-arid Mediterranean. The study area is located in the Mesomediterranean bioclimatic belt [41], which is characterized by a dry climate with scarce rainfalls and a large variability of temperature. The mean annual precipitation is 321 mm and the mean temperature is 16.6 ◦C (minimum and maximum mean temperatures of −2 and 40 ◦C, respectively, in February and July (Spanish Agency for Meteorology, AEMET, period of 1981–2010) [42]. According to the Spanish soil map [43], soils belong to the Aridisol order and Calcic suborder, according to the classification established by the Soil Taxonomy System [44].

A large part of the study area is covered by a sclerophyllous vegetation dominated by *Pinus halepensis* Mill forests and with understory mainly consisting of *Quercus coccifera, Rhamnus lycioides, Halimium hatriplicifolium, Rosmarinus officinalis, Cistus clusii, Rhamnus alaternis, Phamnus alaternis, Genista spartioides* subsp. *Retamoid*; south-facing areas are dominated by *Stipa spartans* and natural therophytic grasslands. Before the wildfire, tree height and density of *Pinus halepensis* M. were 5 to 12 m and 450 to 775 individuals per hectare, respectively [45]. The wildfire occurred on July 2012, and burned a total area of 6500 ha of forest. The fire propagation was very fast, and most of the land was affected by moderate to high fire severity. After the fire, the Forest Service of Castilla-La Mancha region carried out post-fire restoration works, stabilization treatments on hillslopes (both log erosion barriers and contour-felled log debris) in autumn 2012, and checked dams in the catchment reaches in 2013 (Figure 2b,d and Figure 3).

**Figure 2.** Plot location and distribution of soil conditions in the catchment (**a**,**b**) and hillslope (**c**,**d**) areas in Sierra de los Donceles forest (Castilla La Mancha, Spain). Log erosion barriers ("LEB"), contour-felled log debris ("CFD"), burned and no action ("BNA"), unburned ("UB").

**Figure 3.** Photos of the construction of contour-felled log debris (**a**) and log erosion barriers (**b**) in Sierra de los Donceles forest (Castilla La Mancha, Spain).

#### *2.2. Experimental Design*

The catchment selected in the studied area was subject to a wildfire and treated with the two hillslope stabilization treatments, consisting of log erosion barriers ("LEB") and contour-felled log debris ("CFD") (Figure 3). LEBs were built by felling burned trees that are laid on the ground along the slope contour [46]. Each log was anchored in-place and the space between the log and soil surface was filled with soil to create a storage basin upstream of the LEB, where the water and sediment flows are trapped. Earthen berms were sometimes installed to reduce the share of water circumventing the log sides. In the studied catchment, the stabilization treatment was operated at a mean density of 30 LEBs per hectare with a mean length of 10 m (for a linear density of 300 m of logs per ha). These densities were limited by the scarce availability of wood material, due to the unsuitable type of vegetation in the area (small-diameter and low-density trees). The CFD treatment consisted of branch and small-felling burned trees, which were laid on the ground along the slope contour, as for LEB. In this case, logs were not anchored. The mean treatment density was 17 CFD per ha with a mean length of 50 m (corresponding to 850 m per ha), given the lower compaction and concentration of the material for building the CFD.

Moreover, a burned area was left without any treatment (henceforth "burned and no action", BNA) (Figure 2). Another area that was located very close to the burned forest and not affected by fire was considered ("unburned", UB); the vegetation and soil of this area was extremely similar as those of the burned zone and it was therefore representative of the actual pre-fire conditions. Table 1 reports the vegetation cover in the plots under four land conditions at the two spatial scales.


**Table 1.** Pre-fire forest cover in plots under four land conditions before the wildfire of 2012 in Sierra de los Donceles forest (Castilla-La Mancha, Spain).

Notes: UB = unburned; BNA = burned and no action; CFD = contour-felled log debris; LEB = log erosion barriers; data source: surveys of the Forest Service of Castilla-La Mancha region.

The vegetation cover was characterized in all plots under the four soil conditions. More specifically, in the UB plots, the vegetation cover is mainly characterized by woody and herbaceous species (specifically, *Pinus halepensis, Rosmarinus officinalis, Brachypodium retusum* and *Cistus albidus*). *Pinus halepensis, Cistus albidus* and *Rosmarinus officinalis* can largely be disseminated; their seeds, after fire, are stimulated to germinate. Conversely, *Brachypodium retusum* is a herbaceous species (hemicryptophyte), which is a facultative disseminator, with the ability to reproduce both from sprouts and seed and to adapt to frequent fires. In the BNA plots, *Brachypodium retusum, Cistus albidus, Halimium halimifolium, Quercus coccifera, Klasea flavescens* subsp *cichoracea flubensces* ssp. *leucanta* and *Pinus halepensis* are the main species identified after the wildfire. *Quercus coccifera* is a woody phanerophytic species with vegetative propagation. Most of the area in the CFD plots was covered after the wildfire by species with different strategies of fire response, mainly seeding trees (*Cistus albidus, Fumana ericoides, Pinus halepensis*), but also sprouting (*Rhamnus lycioides, Pistacia lentiscus*) and facultative sprouting (*Anthillys cytisoides*) trees as well as facultative sprouting herbaceous species (*Brachypodium retusum*). Finally, also in LEB plots, we surveyed, after the wildfire, a variety of species with different response strategies to the fire. Woody seeding species prevail (*Fumana helicoide, Cistus albidus, Rosmariuns officinalis, Pinus halepensis, Atriplex halimus*), but also woody sprouting trees (*Rhamnus lycioides, Juniperus oxycedrus*), woody facultative sprouting (*Retama sphaerocarpa*), herbaceous facultative sprouting (*Brachypodium retusum, Macrochloa tenacissima* (L.) *Kunth*), and herbaceous seeding (*Asphodelus fistulosus*) species.

#### *2.3. Data Collection and Processing*

#### 2.3.1. Field Survey of Vegetation Cover

Four years after the fire (in the summer of 2016), fifty 30 m × 30 m plots were established in the catchment area (Figure 2a,b) and installed in each of the four land conditions (7 plots in LEB, 7 in CFD, 25 in BNA, and 11 in UB). An additional twelve 10 m × 10 m plots were installed in the hillslope area (Figure 2c,d) for the same land conditions (3 plots in LEB, 3 in CFD, 3 in BNA, and 3 in UB), totaling 62 sample plots (Figure 2c,d). All the plots were separated 100 m from each other, to be considered as independent; moreover, the burned plots were exposed to similar fire severity. In each plot, three strips were selected (10 m long and 0.5 m wide for the 10 m × 10 m plots, and 30 m long and 0.5 m wide for the 30 m × 30 m plots), where the vegetation cover was measured in percentage. Along each strip, we identified the different species, and calculated the percent canopy cover by the line intercept method [46] as the sum of canopy distances by the strip length. The data were averaged in plot as the mean of the three strips [47,48].

#### 2.3.2. Remote Sensing Surveys of Vegetation Cover

On the plots at the hillslope scale, in 2016, a scheduled drone flight was carried out. The UAV used was a quadcopter md4-1000 (Microdrones Inc., Kreuztal, Germany) with a RGB SONY ILCE-5100 digital camera (Sony Corporation, Tokyo, Japan) on board. The sensor of the SONY ILCE-5100 camera was a complementary metal oxide semiconductor (CMOS) Exmor® type APS-C (23.5 × 15.6 mm) with pixel size of 4 × <sup>4</sup> <sup>μ</sup>m. The image size was 6000 × 4000 (columns and rows) and its focal length was 20 mm. Flight planning was performed for a flight altitude of 120 m obtaining a ground sample distance (GSD) of 0.015 m (Figure 4). LANDSAT8 images of the same date were collected over the plots at the catchment scale (Figure 4). *VARI* values were calculated using UAV and LANDSAT8 images, both from 2016. LANDSAT image was subjected to atmospheric correction and radiometric calibration [49] without using LEDAPS approach. *VARI* was calculated for each pixel, which identifies the vegetation in the part of the visible spectrum (ideal for RGB images). *VARI* is calculated as follows [24]:

$$VARI \; = \frac{Green \; - \; Red}{Green \; + \; Red \; - \; Blue} \tag{1}$$

Once calculated, *VARI* was then classified, splitting the entire range of data into five classes, to which a "rank" was given for both image sources (satellite and UAV) (Table 2). Then, possible correlations between the *VARI* values and the vegetation cover measured in the plots in 2016 at the two spatial scales were found. Using the *VARI* values, the effects of the land restoration measures on the vegetation cover in comparison with the burned and not treated, as well as to the unburned land, were evaluated.

**Table 2.** Calculation of the class width for the VARI classification in ranks in Sierra de Los Donceles forest (Castilla-La Mancha, Spain) (Source: our processing).


#### *2.4. Spatial Analysis*

Using QGIS software applied to a digital terrain model (DTM) prepared by the Spanish Center of Geographic Information in 2016 (resolution of 5 × 5 m), a spatial analysis was carried out to calculate the values of land slope and terrain roughness at the hillslope scale. As exposed by Wu et al. [50], terrain roughness is defined as the unevenness of the terrain surface (including rocks and low vegetation) at scales of several meters. This analysis was targeted to identify possible relationships between vegetation regeneration and land characteristics (that is, land slope and terrain roughness) for the different soil conditions. In order to analyze VARI index and its relationships with land slope and terrain roughness at the hillslope scale, a total of 100 plots were randomly selected on the DTM prepared by the Spanish Center of Geographic Information in 2016 (40 plots in LEB, 20 in CFD, 20 in BNA, 20 in UB.

#### *2.5. Statistical Analysis*

For both field measurements and remote sensing estimations, General Linear Models (GLM, using treatment as fixed factor and plot as random factor) was applied to evaluate the statistical significance of the differences in the vegetation cover or VARI among the treatments and the control areas. This statistical approach allows us to deal with pseudoreplication [51]. The homogeneity of the variance and the normality of the samples were checked using the Levene and Kolmogorov-Smirnov tests, respectively. All plots were considered spatially independent. The independent Fisher's least significant difference (LSD) test was used for post hoc comparisons. An α-level <0.05 was adopted. Finally, linear correlations were calculated between VARI and vegetation cover measured in field surveys on one side, as well as land slope and terrain roughness on the other side.

#### **3. Results and Discussion**

#### *3.1. Field Measurements of Vegetation Cover in Different Land Conditions*

All surveyed plant species of this study were typical of post-fire vegetation succession in the Mediterranean forest. These species are commonly found in open areas receiving high solar radiation and adapted to fire through different vegetative and reproductive mechanisms. In more detail, the floristic composition of the shrub and herbal layer in the surveyed area was not significantly altered by fire. In fact, the species that have repopulated the forest areas after fire belonged to the pre-existing populations. The number of species did not change, since after the fire, the species progressively regenerates, thanks to the adaption to new conditions of light, water and nutrients.

The vegetation cover measured in the plots under the four land conditions was more extensive in the areas treated with LEB (vegetal cover of 84.5 ± 4.5%) compared to CFD plots (76.9 ± 5.6%). As expected, the extreme values were surveyed in UB (100%) and BNA (65.5 ± 4.2%) (Figure 5). These differences were significant after one-way ANOVA (*p* < 0.05), except for LEB treatment (not significantly different from UB and CFD). According to fieldwork and forest regional managers, this change is mainly due to shrub regeneration and growth. The Mediterranean vegetation is highly adapted to wildfires and these ecosystems are able to respond to fire using resprout or seeding mechanisms. That change is attributed to the higher biomass of the scrub and herbal layer. As recently stated by Keeley and Pausas [52], fire does not threaten ecosystem health, since it is a necessary process and a natural disturbance that is beneficial for the functionality of a fire-adapted ecosystem in Mediterranean forests. As demonstrated by our results, following fire disturbance, Mediterranean forests can regenerate from seed banks stored in soils and basal resprouts.

**Figure 5.** Vegetation cover (mean ± standard deviation of 62 plots) in plots under four land conditions (UB = unburned; BNA = burned and no action; CFD = contour-felled log debris; LEB = log erosion barriers) after the wildfire of 2012 in Sierra de Los Donceles forest (Castilla-La Mancha, Spain). Mean values that do not share a lower case letter (top of graph) are significantly different from each other (HSD, *p* < 0.05).

This means that, four years after the wildfire, the vegetation regeneration is far from covering the entire area, such as in the unburned zone. The implementation of post-fire treatments allows for a significantly faster recovery process, and LEBs are particularly effective. This higher regeneration compared to CFD-treated areas is somewhat expected, since the dead and burned material of which CFDs consist of is degrading, as well as incorporated more easily compared to the burned wood of LEBs [53]. Both post-fire management techniques are recognized to reduce runoff and erosion (by slowing flows of water and sediments) [14] and improve the physico-chemical properties of soils (increasing infiltration and water retention, and organic matter and nutrients), thus enhancing establishment and development of vegetation [12]. Some LEBs can trap up to 40% of sediments; moreover, this treatment is cheaper compared to other hillslope stabilization techniques [54]. CFDs are effective to reduce water and sediment flows in burned forest subject to machinery salvage logging [55]. Moreover, according to Badía et al. [13], hillslope stabilization after wildfire is a physical barrier that avoids losses or reductions in soil organic matter and nutrients. At the same time, logs of CFDs and canopy residues of LEBs can change the microclimate and land conditions, as well as represent organic matter source after decomposition. This enhances the biological activity of soils [18]. Moreover, decreased evaporation, higher soil moisture, and soil organic matter accumulation upslope of LEB and CFD might lead to an increase in soil respiration and microbial activity, also enhancing nutrient availability after a burn at 5 years from the wildfire). Regarding other studies evaluating post-fire regeneration of vegetation, Christakopoulos et al. [40] found that reforestation was, in some cases, comparable, and in other areas, higher, compared to natural regeneration of pines in burned forests of Greece.

#### *3.2. Correlations of VARI with Vegetation Cover Using Remote Sensing Techniques*

Figure 6a,b reports the maps of VARI distribution in the experimental areas at both catchment and hillslope scales. In particular, Figure 6b highlights that VARI is higher in the areas subjected to post-fire treatments and lower in burned and untreated land.

**Figure 6.** Spatial distribution of VARI surveyed by satellite (**a**) and UAV (**b**) images of 2016 among four land conditions after the wildfire of 2012 in Sierra de Los Donceles forest (Castilla-La Mancha, Spain). Legend: Fajinas = log erosion barriers; cordones = contour-felled log debris; sin tratamiento = burned and no Action; sin fuego = unburned.

A low and not significant R2 (0.03, *p* > 0.05) was found regressing the LANDSAT8 derived VARI with the vegetation cover measured in field (Figure 7a, catchment scale), while the regression of UAV-derived VARI versus vegetation cover was significant (R2 = 0.84, *p* < 0.05) (Figure 7b, hillslope scale). This difference is clearly due to the lower spatial resolution of satellite data compared to surveys by UAV.

**Figure 7.** Correlations between VARI and vegetation cover at catchment (**a**), and hillslope (**b**) scales surveyed by LANDSAT8 and UAV images (2016) in Sierra de Los Donceles forest (Castilla-La Mancha, Spain).

The regression equations for the correlation VARI versus vegetal cover (VC) measured in field was the following:

$$\text{VC} = 16.90 \text{ VARI} + 9.15 \tag{2}$$

The vegetation cover distribution in classes according to the VARI classification by UAV is reported in Table 2 and confirms the spatial differences among the land conditions. By visually comparing the vegetation cover distribution according to the two VARI classifications (satellite and UAV, Table 2), no overlay among the same class was detected. This lack of correspondence further confirms the higher reliability of VARI estimated from UAVs to reproduce land cover compared to estimations by satellite, due to unsuitability of image resolution.

The use of VARI to evaluate vegetation regeneration from UAVs has been widely studied in literature. To cite the most recent studies, VARI was used to analyze vegetation on different land uses (e.g., [36–39,56]), monitoring vegetation into water bodies (e.g., [57], and preparation of the Digital Elevation Model (DEM) [58]). This index shows significant correlation with crop height and yield [59,60]. In burned lands, Larrinaga and Brotons [39] calculated VARI for analysis of post fire regeneration of Mediterranean forests. However, these authors reported that this index underperformed compared to other greenness indexes, such as the Excess green index (ExGI) and green chromatic coordinate (GCC) index, and the same was found for the green red vegetation index (GRVI). Despite this, these authors have demonstrated that low-cost UAVs may improve forest monitoring after disturbance, even in those habitats and situations where resource limitation is an issue. In general, VARI shows a minimal sensitivity to atmospheric effects [26,61,62].

#### *3.3. Evaluation of the Vegetation Regeneration in Fire-Affected and Treated Areas Using VARI*

The comparison of the mean VARI values under the four studied land conditions between the two remote images shows that, at the catchment scale, VARI was higher (0.009 ± 0.046) under CFDs and lower (−0.034 ± 0.035) in UB areas. This appears again unrealistic, further confirming the misleading meaning of land cover images captured by UAV. Moreover, the differences among the four land conditions were not significant (*p* = 0.724) (Figure 8a).

**Figure 8.** Mean values of VARI surveyed by LANDSAT8 (**a**) and UAV (**b**) images (2016) among four land conditions after the wildfire of 2012 in Sierra de Los Donceles forest (Castilla-La Mancha, Spain).

In contrast, when analyzing the images by UAV, the highest VARI (4.578 ± 0.230) was found in the UB areas (indicating the most extensive vegetation cover), while the lowest value was detected in the areas burned but not treated (VARI of 3.024 ± 0.853). However, in contrast with field surveys, vegetation regeneration was higher in the land treated with CFDs (mean VARI of 4.270 ± 0.350) compared to the areas with LEBs (3.889 ± 0.416) (Figure 8b). As for field surveys, the differences in VARI were significant (*p* < 0.05) at this smaller scale using the UAV. Lack of correspondence between the post-fire land treatments suggests that UAV is a viable technique in detecting the variability of vegetation cover only when the difference in VARI mean values is noticeable (e.g., between burned and unburned areas); conversely, the contrasts between areas with similar VARIs may be quite misleading, for instance, when the effectiveness of a treatment has to be evaluated. This result is in close accordance with the indications by Corona et al. [63], who suggested the use of aerial and satellite imagery characterized by high or, better, very high spatial resolution for an effective support to post-fire management (burned area mapping, fire severity assessment, post-fire vegetation monitoring).

Comparisons to other published data about evaluations of post-fire vegetation recovery using remote sensing and greenness indices reveal that SAR images perform better compared to optical images in estimating forest regeneration, particularly when objectbased classification procedure is applied (nearly 90% of accuracy) [34]. Moreover, Lasaponara et al. [33] found a better reliability of NDVI (Normalized Difference Vegetation Index) in capturing the diverse vegetation regeneration in both natural and managed areas as well as before and after fire occurrence compared to NBR index, thanks to the data processing using detrended fluctuation analysis (DFA). LANDSAT time series analysis with NBR application was a useful means of describing and analyzing post-fire vegetation recovery across mixed-severity wildfire extents [30], while, in contrast, Lentile et al. [31] reported that dNBR is an imperfect indicator of post-fire effects on vegetation and soil, when the burned areas are affected by different burn severities and are of a highly variable patch size. Regarding the use of drone images, Fernández-Guisuraga et al. [64] highlighted the more detailed spatial information provided by the drone orthomosaic compared to WorldView-2 satellite imagery in vegetation regeneration of heterogeneous burned areas.

Overall, the use of UAV systems for vegetation recovery monitoring is still under development, since these systems have a great potential not yet fully valorized; however,

this use is still limited by some constraints, such as battery life and availability of cameras with suitable spectral range, as well as cost [22].

#### *3.4. Spatial Distribution of Land Slope and Roughness, and Correlations with VARI*

The spatial analysis carried out in the experimental areas showed that the slope of land subjected to post-fire treatments was significantly different (32.2 ± 5.26% for LEB and 36.2 ± 2.61% for CFD) compared to the BNA (23.1 ± 6.22%), but similar to the UB land (36.6 ± 3.71%) (Table 3 and Figure 9a). This may be due to the fact that landscape managers prioritized steeper areas when post-fire treatments were planned, since the zones with higher slopes are more prone to erosion and hydrogeological risks compared to flat areas.

**Table 3.** Values (mean ± standard deviation) of land slope and terrain roughness in the small-scale areas under four land conditions in Sierra de los Donceles forest (Castilla-La Mancha, Spain).


Notes: UB = unburned; BNA = burned and no action; CFD = contour-felled log debris; LEB = log erosion barriers; different lower case letters indicate significant differences (HSD, *p* < 0.05).

**Figure 9.** *Cont*.

**Figure 9.** Spatial distribution of land slope (**a**) and terrain roughness (**b**) surveyed by UAV images of 2016 among four land conditions after the wildfire of 2012 in Sierra de Los Donceles forest (Castilla-La Mancha, Spain). Legend: Fajinas = log erosion barriers; cordones = contour-felled log debris; sin tratamiento = burned and no action; sin fuego = unburned.

According to the land map of Figure 9b prepared using spatial analysis, the terrain roughness of the burned areas was significantly different (0.08 ± 0.03 μm for LEB, 0.09 ± 0.01 μm for CFD, and 0.06 ± 0.03 μm for BNA) compared to UB zone (0.28 ± 0.21 μm) (Table 2). This means that the vegetation growing in the burned and treated zone do not have a significant effect on the terrain roughness compared to the area without treatment, and that these values are well below the roughness of the unburned soil. A high terrain roughness is beneficial to reduce surface runoff and thus, sediment transport downstream [65–67]. Soil preparation (e.g., by tillage, conditioning, and terracing) after wildfire may be suggested, in order to enhance vegetation growth and improve the soil's resistance to erosion [68–71].

A fair and significant linear correlation was found regressing VARI on land slope under all soil conditions (r<sup>2</sup> = 0.52, *p* < 0.05) (Figure 10a), while the linear regression between VARI and terrain roughness was much lower (r<sup>2</sup> = 0.15) and not significant (*p* = 0.49) (Figure 10b). This correlation becomes significant (*p* < 0.05) although, again, low (r2 = 0.35 to 0.38), adopting exponential or logarithmic equations, whose physical meaning is, however, difficult to be justified (data not shown).

**Figure 10.** Scatterplots of vegetation regeneration (measured by VARI by UAV images of 2016) versus land slope (**a**) and terrain roughness (**b**) among four land conditions after the wildfire of 2012 in Sierra de Los Donceles forest (Castilla-La Mancha, Spain).

The first correlation suggests that the effectiveness of post-fire treatments on vegetation regeneration (measured using VARI) increased with slope, and thus, the strategy of landscape managers that have prioritized steeper slopes for treatments was successful. By contrast, the lower, or absence of, significance of correlations between VARI and terrain roughness implies that vegetation growth, at least in the experimental conditions, played a limited influence on soil surface conditions. However, vegetation regeneration remains a key factor to reduce soil exposure to erosion in wildfire-affected areas.

#### **4. Conclusions**

The surveys of vegetation cover using fieldwork and remote sensing in lands burned by wildfire in a Mediterranean forest showed that:


The spatial analysis of distribution of VARI and land characteristics at the hillslope scale proved that:


However, the validity of these preliminary results must be confirmed in other experimental conditions, such as in areas with different climate and soil conditions, as well as in soils subjected to other post-fire management techniques. The effects of patch vegetation should be also explored, since a high spatial variability of vegetation may not be fully captured by aerial or satellite images. Additionally, the implementation of UAV-based surveys throughout the different stages of vegetation regeneration after wildfire can be suggested, in order to develop a viable tool for monitoring the effectiveness of post-fire techniques over time. A wider validation activity should ensure a practical use of UAV images to support the activity of land managers in planning and implementing efficient measures to restore wildfire-affected areas using vegetation cover. The research question about the viability of higher-resolution satellite images to detect contrasts in vegetation cover among different land conditions remains open. Presumably, the integration of several techniques, combining advantages and limiting constraints, can be suggested for estimation of post-fire vegetation recovery in forest ecosystems with dynamically variable characteristics.

**Author Contributions:** Conceptualization, M.E.L.-B. and P.A.P.-A.; methodology, M.E.L.-B., M.A.M., D.H., J.G.-R., D.A.Z., and P.A.P.-A.; formal analysis, J.L.M. and P.D.; writing—original draft preparation, J.L.M., M.E.L.-B., P.A.P.-A., P.D., M.A.M., D.H., J.G.-R., and D.A.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Superpixel-Based Singular Spectrum Analysis for Effective Spatial-Spectral Feature Extraction**

**Subhashree Subudhi 1, Ramnarayan Patro 1, Pradyut Kumar Biswal 1,\* and Fabio Dell'Acqua <sup>2</sup>**


**Abstract:** In the processing of remotely sensed data, classification may be preceded by feature extraction, which helps in making the most informative parts of the data emerge. Effective feature extraction may boost the efficiency and accuracy of the following classification, and hence various methods have been proposed to perform it. Recently, Singular Spectrum Analysis (SSA) and its 2-D variation (2D-SSA) have emerged as popular, cutting-edge technologies for effective feature extraction in Hyperspectral Images (HSI). Using 2D-SSA, each band image of an HSI is initially decomposed into various components, and then the image is reconstructed using the most significant eigen-tuples relative to their eigen-values, which represent strong spatial features for the classification task. However, instead of performing reconstruction on the whole image, it may be more effective to apply reconstruction to object-specific spatial regions, which is the proposed objective of this research. As an HSI may cover a large area, multiple objects are generally present within a single scene. Hence, spatial information can be highlighted accurately by specializing the reconstruction based on the local context. The local context may be defined by the so-called superpixels, i.e., finite sets of pixels that constitute a homogeneous set. Each superpixel may undergo tailored reconstruction, with a process expected to perform better than non-spatially-adaptive approaches. In this paper, a Superpixel-based SSA (SP-SSA) method is proposed where the image is first segmented into multiple regions using a superpixel segmentation approach. Next, each segment is individually reconstructed using 2D-SSA. In doing so, the spatial contextual information is preserved, leading to better classifier performance. The performance of the reconstructed features is evaluated using an SVM classifier. Experiments on four popular benchmark datasets reveal that, in terms of the classification accuracy, the proposed approach overperforms the standard SSA technique and various common spatio-spectral classification methods.

**Keywords:** hyperspectral image; superpixel segmentation; evaluation; 2D-singular spectrum analysis (2D-SSA); feature extraction

#### **1. Introduction**

Recent advancements in hyperspectral sensors resulted in the increased availability of Hyperspectral Images (HSI) and a boost in their circulation among the remote sensing community. HSI data enables the discrimination of objects even with minor differences as it contains several contiguous spectral bands acquired from the visible to the infrared region [1] so that every small spectral difference can, in principle, be captured. The information is available in the form of a 3-D structure that contains a 2-D spatial scene along with a 1-D spectral signature. These unique characteristics of HSI have made them popular in several application areas, such as agriculture [2], mineralogy [3], land cover classification [4], target detection [5], and others. However, effective classification of HSI is still an open challenge.

**Citation:** Subudhi, S.; Patro, R.N.; Biswal, P.K.; Dell'Acqua, F. Superpixel-Based Singular Spectrum Analysis for Effective Spatial-Spectral Feature Extraction. *Appl. Sci.* **2021**, *11*, 10876. https://doi.org/10.3390/ app112210876

Academic Editors: Dario Gioia and Maria Danese

Received: 29 September 2021 Accepted: 2 November 2021 Published: 17 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Several classification techniques, such as K-nearest neighbor (KNN) [6], support vector machine (SVM) [7], multinomial logistic regression (MLR) [8], Extreme Learning Machine (ELM) [9], and Sparse Representation Classifier (SRC) [10] have been proposed in the past decades. The richness in spectral information attracted research efforts on pixel-based processing and classification. SVM is the most popular and widespread classifier due to its lower generalization error rate that makes it capable of identifying even minor changes in spectral signatures. Due to the high spectral dimensionality compared to a generally limited number of class-specific training samples, it is quite difficult to properly estimate the model parameters. Hence, there is a need to adopt effective spatial-spectral feature extraction approaches to overcome the aforementioned challenges.

To deal with the issue of higher spectral dimensionality, several linear (e.g., Principal Component Analysis (PCA) [11], Independent Component Analysis (ICA) [12], and Linear Discriminant Analysis (LDA) [13]) and non-linear (manifold learning [14]) dimensionality reduction (DR) methods have been introduced. Band selection approaches may also be utilized to select the most informative bands out of several available bands [15,16]. Spectral features alone, however, may not be sufficient to score very high accuracy values. To improve performance, it is necessary to incorporate spatio-spectral features that help increasing the separability of classes.

Recently, various spatial feature extraction techniques have been proposed. Mathematical Morphology is one of the most popular approaches that is extensively utilized by researchers. The concept of Extended Morphological Profiles (EMP) for FE in HSI was first proposed by Benediktsson et al. [17]. This technique utilizes morphological opening and closing transformations to extract spatial geometrical information. Later, Dalla Mura et al. [18] proposed Morphological Attribute Filters (MAP) for the spatial FE. From that point onwards, several variations of Attribute Profiles (AP) were created. Ghamisi et al. [19] conducted a comprehensive survey on the evolution in Attribute Profiles.

Texture Descriptors, including Wavelet transform [20], Gray-Level Co-occurrence Matrix (GLCM) [21], Local Binary Patterns (LBP) [22], and Gabor filters [23], have also been used in the literature for spatial FE. Filtering, i.e., moving-window-based processing, is another approach to extract spatial-spectral features. Various edge-preserving filters, such as Bilateral Filters [24], Trilateral Filters [25], Guidance filters [26], and Domain Transform Recursive Filters [27] have been introduced for spatial FE in the literature. The texture and noise variations are minimized by performing smoothing operations; however, important details, such as edges and lines, are well preserved by these filters [28].

In addition to these techniques, 2D-SSA is another interesting approach for spatial feature extraction. Using this approach, each band image of HSI is initially decomposed into varying trends, oscillations, and noise. Later, HSI is reconstructed using the selected oscillations and trends [29,30]. In 2D-SSA, the spatial structural information is extracted by utilizing the characteristics of the surrounding pixels in a specific embedding window. It can withstand high levels of noise and generally achieves good data classification results.

2D-SSA suffers, however, from several limitations, such as reduced utilization of the abundant spectral information available in data, and over-smoothing or under-smoothing of classification results because of the fixed embedding window size. To overcome the challenge of selecting the optimal embedding window size, recently, multi-scale 2D-SSA has been proposed for the effective extraction of discriminative spatial features under different noise conditions [31,32].

Superpixel segmentation techniques have gained popularity in recent years due to their capability of exploiting spatial structural information adaptively in an image. In [33], a survey on superpixel segmentation as a preprocessing step in HSI analysis was presented. A superpixel-based classification via multiple kernels (SCMK) approach was proposed in [34]. In [35], a region-based relaxed multiple kernel (R2MK) method was proposed that combines the multiscale spectral and spatial features using a kernel collaborative representation classification technique.

To obtain superior classification performance and solve the problem of optimal superpixel number selection, an adjacent superpixel-based multiscale spatial-spectral kernel (ASMGSSK) was proposed in [36]. In [37], a multiscale segmentation-based SuperPCA model (MSuperPCA) was developed, which can effectively integrate multiscale spatial information to obtain the optimal classification result by decision fusion.

Recently, deep learning techniques have become quite popular in the classification of HSI data due to their ability to extract discriminant and abstract features by using a series of hierarchical layers. The initial layers usually extract texture and edge information, whereas deeper layers highlight more complicated features. Some of the most popular deep learning frameworks include stacked autoencoders (SAE) [38], Deep Belief Networks (DBN) [39], Convolution Neural Networks (CNN) [40], Recurrent Neural Networks (RNN) [41], Generative Adversarial Networks (GAN) [42], etc.

Although deep learning approaches have several advantages, they also pose significant challenges in HSI applications. First of all, to achieve better classification result, often deep learning techniques demand large volumes of training samples. Moreover, a large number of hyper-parameters (like the kernel sizes, learning rate, etc.) are involved in training complex deep learning networks mainly designed for feature extraction and classification. Hence, the process becomes computationally expensive.

The disadvantages of combining SSA with structured approaches to incorporating spatial information may be overcome by using more flexible ways to spatially partition the dataset. In line with this consideration, in this work, a superpixel-based SSA (SP-SSA) algorithm was proposed as a means to increase the classifier performance. Instead of performing direct reconstruction, an object-specific reconstruction is performed to accurately preserve the local contextual information. Superpixel segmentation is first applied on the input HSI to generate a segmented HSI where each sub-region carries similar characteristic features, and its shape and size is adjusted according to the local image structure information. Next, 2D-SSA is individually applied on each segmented region to produce the reconstructed HSI. Lastly, the final classification map is generated by using the popular SVM classifier. The major novel contributions of this work are highlighted in the following list:


The remainder of this paper is organized as follows. A detailed description of the proposed method is presented in Section 2. The experimental setup, results, and analysis are described in Section 3. Finally, some conclusions and future work are discussed in Section 4.

#### **2. The Proposed Methodology**

The proposed SP-SSA method includes three stages as described in Figure 1. In stage 1, superpixel segmentation is applied on the input HSI to obtain the segmented HSI. In stage 2, each segmented region is reconstructed using 2D-SSA to obtain the reconstructed HSI. In the final stage, an SVM classifier is applied on the reconstructed HSI to build the final classification map. A detailed description of each of these stages is presented in the subsections below.

**Figure 1.** Flowchart of the proposed method.

#### *2.1. Superpixel Segmentation*

Superpixel segmentation approaches have gained popularity in recent years as these approaches have several benefits. Using superpixels, the computational complexity can be drastically reduced by computing features on more meaningful regions rather than acting on each individual pixel in HSI [43]. Simple Linear Iterative Clustering (SLIC) [44] is one of the most popular gradient-ascent-based superpixel segmentation approaches, where an initially defined tentative set of cluster points are iteratively refined using a gradient-ascent method until some convergence criteria are met. This algorithm has lower computational complexity as it applies the *k*-means method locally. The algorithm includes four key steps that can be summarized as follows.

The first step is cluster center initialization. Let the input HSI be denoted as *<sup>H</sup><sup>b</sup>* ≡ {*hb* <sup>1</sup>, *<sup>h</sup><sup>b</sup>* <sup>2</sup>, ... , *<sup>h</sup><sup>b</sup> <sup>N</sup>*} with *<sup>N</sup>* pixels, where {*h<sup>b</sup> <sup>i</sup>* } represents the value at the *i*th pixel for the *b*th spectral band and *i* = 1, 2, ... *N*; *b* = 1, 2, ... *B*. *B* is the total number of spectral bands. Each pixel can be labeled as *Ai* = [*hi*,*ri*, *ui*], where *h<sup>T</sup> <sup>i</sup>* = [*h*1, *h*2, ... , *hB*] *<sup>T</sup>* is the spectral vector and [*ri*, *ui*] *<sup>T</sup>* is the position vector. The *<sup>K</sup>* number of initial cluster centers C*<sup>j</sup>* = [*hj*,*rj*, *uj*] *T* -

are sampled on a regular *Q* × *Q* (*Q* = *N <sup>K</sup>* ) grid and are, thus, equally spaced apart [45].

The next step is the cluster assignment step, where each pixel is assigned to the nearby cluster center based on the computed distance measure *D*. Distance is computed within a 2*Q* × 2*Q* window around the cluster center. The distance between the cluster center C*<sup>j</sup>* and pixel *Ai* is calculated as follows (Equation (1)):

$$D = D\_{\text{spectral}} + \frac{w}{Q} D\_{\text{spatial}} \tag{1}$$

where *w* is the weighting factor between spectral and spatial features. The spectral and spatial distance between pixel *i* and *j* are represented as in Equations (2) and (3) below.

$$D\_{\text{spectral}} = \sqrt{\sum\_{b=1}^{B} \left( h\_i^b - h\_j^b \right)^2} \tag{2}$$

where *Dspectral* is the measure of homogeneity within the superpixels.

$$D\_{\text{spatial}} = \sqrt{(r\_i - r\_j)^2 + (u\_i - u\_j)^2} \tag{3}$$

where (*ri*, *ui*) denotes the location of pixel *i* in the superpixels. The spatial distance *Dspatial* ensures regularity and compactness in the generated superpixels.

In the third step, the cluster centers are updated with the mean value of all pixels belonging to the same cluster. The second and third steps are iteratively repeated until convergence is achieved.

In the final step, post-processing is performed to enforce connectivity by reassigning disjoint pixels to nearby superpixels.

#### *2.2. 2D-SSA*

SSA is capable of decomposing a series into multiple independent components or subseries, where each extracted eigenvalue represents an individual component of the

original series. The SSA can be applied to the respective spectral bands of the hypercube, thereby, decomposing the 2-D scene, and then reconstructing it using the respective main components while removing the noise contribution. As a data cube is decomposed in this way, the local structure and main spatial trends are typically found in the first component. Hence, when all images within the hyperspectral cube are decomposed and only the first components are selected to individually reconstruct each of them, a resulting cube with minimum noise is generated. The SSA can be implemented using the following four steps:

#### 2.2.1. Embedding

Imagine a HSI dataset *H*, with a size of *Nx* × *Ny* × *B*, where *Nx*, *Ny* indicates the band image size and *B* represents the total number of available bands. Each band image *H<sup>b</sup>* (*b* ∈ *B*) can be expressed as follows:

$$H^b = \begin{bmatrix} H\_{1,1}^b & H\_{1,2}^b & \cdots & H\_{1,N\_y}^b \\ H\_{2,1}^b & H\_{2,2}^b & \cdots & H\_{2,N\_y}^b \\ \vdots & \vdots & \ddots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ H\_{N\_x,1}^b & H\_{N\_x,2}^b & \cdots & H\_{N\_x,N\_y}^b \end{bmatrix}\_{N\_x \times N\_y} \tag{4}$$

Next, a 2D window *<sup>Q</sup><sup>b</sup>* is defined, whose dimensions are *Mx* × *My*.

$$Q^b = \begin{bmatrix} H\_{i,j}^b & H\_{i,j+1}^b & \cdots & H\_{i,j+M\_y-1}^b \\ H\_{i+1,j}^b & H\_{i+1,j+1}^b & \cdots & H\_{i+1,j+M\_y-1}^b \\ \vdots & \vdots & \ddots & \vdots \\ H\_{i+M\_x-1,j}^b H\_{i+M\_x-1,j+1}^b & \cdots & H\_{i+M\_x-1,j+M\_y-1}^b \end{bmatrix}\_{M\_x \times M\_y} \tag{5}$$

where 1 ≤ *Mx* ≤ *Nx*, 1 ≤ *My* ≤ *Ny*, and 1 < *MxMy* < *NxNy*. Each pixel is spatially positioned by (*i*, *j*) within the image *Hb*. The pixels in a window *Q<sup>b</sup>* can be rearranged into a column vector *C<sup>b</sup> <sup>i</sup>*,*<sup>j</sup>* <sup>∈</sup> <sup>R</sup>*MxMy* according to the reference position (*i*, *<sup>j</sup>*) as follows:

$$C\_{i,j}^b = [H\_{i,j'}^b H\_{i,j+1'}^b, \dots, H\_{i,j+M\_Y-1'}^b, H\_{i+1,j'}^b, H\_{i+1,j+1'}^b, \dots, H\_{i+M\_X-1,j+M\_Y-1}^b]^T \tag{6}$$

To scan the whole image *Hb*, this 2-D window is slid across it from top left to bottom right until it has visited every position on the entire image (see also Figure 2 for a graphical explanation).

**Figure 2.** Moving window across the image *H<sup>b</sup>* to create the trajectory matrix Z*b*.

As a result, the trajectory matrix Z*<sup>b</sup>* of all feasible 2-D windows of image *H<sup>b</sup>* of size *MxMy* × (*Nx* − *Mx* + 1)(*Ny* − *My* + 1) can be obtained as follows:

$$\mathbf{G}^{b} = \left[ (\mathbf{C}\_{1,1}^{b})^{T}, (\mathbf{C}\_{1,2}^{b})^{T}, \dots, (\mathbf{C}\_{1,N\_{\mathcal{Y}}-M\_{\mathcal{Y}}+1}^{b})^{T}, (\mathbf{C}\_{2,1}^{b})^{T}, \dots, (\mathbf{C}\_{N\_{\mathcal{X}}-M\_{\mathcal{X}}+1, N\_{\mathcal{Y}}-M\_{\mathcal{Y}}+1}^{b})^{T} \right]\_{M\_{\mathcal{X}}M\_{\mathcal{Y}} \times (N\_{\mathcal{X}}-M\_{\mathcal{X}}+1)(N\_{\mathcal{Y}}-M\_{\mathcal{Y}}+1)} \tag{7}$$

Note that the trajectory matrix Z*<sup>b</sup>* has a structure of Hankel–block–Hankel (HbH). Z*<sup>b</sup>* can be expressed as follows:

$$\mathbf{G}^{b} = \begin{bmatrix} P\_1^b & P\_2^b & \cdots & P\_{N\_x - M\_x + 1}^b \\ P\_2^b & P\_3^b & \cdots & P\_{N\_x - M\_x + 2}^b \\ \vdots & \vdots & \ddots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ P\_{M\_x}^b & P\_{M\_x + 1}^b & \cdots & P\_{N\_x}^b \end{bmatrix} \tag{8}$$

Each of the submatrices *P<sup>b</sup> <sup>i</sup>* corresponds to a Hankel structure as follows:

$$P\_i^b = \begin{bmatrix} H\_{i,1}^b & H\_{i,2}^b & \cdots & H\_{i,N\_y-M\_Y+1}^b \\ H\_{i,2}^b & H\_{i,3}^b & \cdots & H\_{i,N\_y-M\_Y+1}^b \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ H\_{i,M\_y}^b & H\_{i,M\_y+1}^b & \cdot & \cdot \end{bmatrix}\_{M\_y \times (N\_y-M\_Y+1)}$$

#### 2.2.2. Singular Value Decomposition (SVD)

After obtaining the trajectory matrix Z*b*, SVD is applied to determine the eigenvalues *λ*<sup>1</sup> ≥ *λ*<sup>2</sup> ≥···≥ *λMxMy* , and the corresponding eigenvectors *U*1, *U*2, ···, *UMxMy* of Z*b* Z*b T* . It is possible to rewrite Z*<sup>b</sup>* as follows:

$$\mathfrak{Z}^{\rho} = \mathfrak{Z}\_1^{\rho} + \mathfrak{Z}\_2^{\rho} + \dots + \mathfrak{Z}\_{M\_x M\_y}^{\rho} \tag{10}$$

where the *i*th elementary matrix is Z*<sup>b</sup> <sup>i</sup>* <sup>=</sup> <sup>√</sup>*λiUiV<sup>T</sup> <sup>i</sup>* and its Principal Components (PCs) are defined as *Vi* <sup>=</sup> (Z*b*) <sup>√</sup> *TUi λi*

#### 2.2.3. Grouping

A subsequent operation is eigenvalue grouping, during which the total set of *MxMy* individual components in (10) are divided into *m* subsets, designated as *S* = [*S*1, *S*2, ... , *Sm*]. By selecting one or more elementary matrices Z*<sup>b</sup> <sup>i</sup>* from each subset, it is possible to derive the main information contained in an image without being disturbed by high noise levels. As a result, the trajectory matrix Z*<sup>b</sup>* can be represented as follows:

$$\mathfrak{Z}^b = \mathfrak{Z}\_{\mathbb{S}\_1}^b + \mathfrak{Z}\_{\mathbb{S}\_2}^b + \dots + \mathfrak{Z}\_{\mathbb{S}\_m}^b \tag{11}$$

The reconstruction of a single band scene of HSI using various numbers of components (Z<sup>b</sup> <sup>i</sup> ) is compared in Figure 3. In general, the component with the highest eigenvalue is the most informative one, containing key features with the lowest noise contribution. With the inclusion of additional components, the reconstructed scene begins to resemble the actual scene. The reconstructed image obtained by grouping the 1st–5th components and 1st–10th components are very similar with marginal differences (Figure 3c,d). Hence, a small number of key components are sufficient to reconstruct the scene satisfactorily.

**Figure 3.** Implementation of 2D-SSA on a HSI scene (**a**) Original scene at 667 nm. (**b**) 1st component grouping. (**c**) 1–5th component grouping. (**d**) 1–10th component grouping, where *Mx* = 5, and *My* = 5.

#### 2.2.4. Diagonal Averaging

Z*b*, in this case, does not necessarily belong to the HbH matrix type. It is projected into a 2D-signal by applying the Hankelization process in two steps; first inside every block (9) and next block-to-block (8) by averaging the anti-diagonal elements in the matrix. Thus, it is possible to obtain a reconstructed image that contains the distinctive spatial features based on the local contextual information present in a 2D window defined by the user.

#### *2.3. Novelty of the Proposed SP-SSA Method*

The proposed approach integrates SSA and superpixel segmentation for the first time to extract improved the spatio-spectral features from HSI. Reconstruction of objectspecific spatial sections, rather than the entire image, may be more effective. Hence, in the proposed work, 2D-SSA is applied individually to each superpixel segmented region to extract the local contextual information accurately. The pseudo-code for the proposed SP-SSA algorithm is outlined in Algorithm 1.


For each superpixel, the reconstruction (*reconstruct*2*DSSA* (Algorithm 1)) is applied to the rectangular Region of Interest (ROI) surrounding the superpixel (Figure 4). The ROI is created based on the location information of the pixels available in that particular segment. Only the reconstructed pixels specific to those pixels in the selected superpixel are stored as spatial features, while the remaining reconstructed pixels in the ROI are discarded as they do not belong to the superpixel under test. The same procedure is applied to all other superpixels, and the HS image is reconstructed using the proposed SP-SSA approach. This procedure collects local object-specific superpixel-based spatial features for each band in the image.

#### *2.4. Classification*

The selection of an appropriate classifier is critical in assessing the performance of the above-mentioned features, especially in hyperspectral images with a limited number of training samples. SVM is the most widely used supervised statistical learning framework among pixel-wise classifiers. With the help of a kernel function, data can be mapped to a higher-dimensional space via a nonlinear transformation, aiming to determine the best hyperplane for separating samples belonging to different classes. The performance of SVM in HSI classification is outstanding despite the variation of the data dimensions [46,47]. Hence, in this work, the SVM classifier is utilized to evaluate the performance of the reconstructed features.

#### **3. Results and Discussion**

This section reports the outcome of testing the proposed approach on some of the most popular benchmark datasets and compares it with other, state-of-art classification approaches.

#### *3.1. Dataset Description*

In this subsection, the datasets used for testing the proposed approach are presented and described.

#### 3.1.1. Indian Pines

The first dataset, named "Indian Pines" (IP), was collected over Northwestern Indiana, USA, with the airborne AVIRIS sensor; it includes a total of 220 bands covering wavelengths from 0.4 to 2.5 μm. About 70% of the imaged area is agricultural land, while the remaining portions are forests. Due to the comparatively low spatial resolution (20 m/pixel) of the sensor, this dataset is challenging as it contains highly mixed pixels. The number of samples obtained per class is also unbalanced, which further complicates classification. The size of the scene is 145 × 145 pixels, and its Ground Truth (GT) data defines 16 different classes. The pseudo-color image, the GT map, and the class names for the dataset are all included in Figure 5.

**Figure 5.** (**a**) False Color Composite Image, (**b**) Ground Truth Image and (**c**) Class names for the Indian Pines Dataset.

#### 3.1.2. Pavia University

The ROSIS sensor was instrumental to the collection of this dataset over the University of Pavia, Italy. The dataset is called "Pavia University" (PU). It has a spatial resolution of 1.3 m and originally comprises 115 spectral bands covering wavelength ranges from 0.43 to 0.86 μm. In the final analysis, 103 bands are used after the elimination of noisy channels. The image has a size of 610 × 340 pixels, and it has nine challenging classes with nearly similar spectral reflectances. Detailed information about the false-color image, Ground Truth, and class names is displayed in Figure 6.

**Figure 6.** (**a**) False Color Composite Image, (**b**) Ground Truth Image and (**c**) Class names for the Pavia University Dataset.

#### 3.1.3. Salinas Dataset

The "Salinas" (SAL) dataset was captured over the Salinas Valley, California, USA, using the AVIRIS Sensor. The sensor has 224 channels with spectral range varying from 0.43 μm to 2.5 μm. This scene has a size of 512 ×217 pixels and spatial resolution of 3.7 m per pixel. The number of bands reduces to 204 after discarding 20 water absorption bands: [108–112], [154–167], 224. The scene is mainly an agricultural area, with 16 classes in its Ground Truth. A false color representation, the Ground Truth, and the class names for the Salinas dataset are shown in Figure 7.

**Figure 7.** (**a**) False Color Composite Image, (**b**) Ground Truth Image, and (**c**) Class names for the Salinas Dataset.

#### 3.1.4. Houston 2018

The 2018 IEEE GRSS Data Fusion Contest (DFC) triggered public dissemination of this rich dataset, which was included in our tests to increase their statistical significance. The image of the Houston campus and its surrounding area was captured by the IRTES CASI-1500 sensor at a GSD of 1 m over Houston, Texas, USA. It has 601 × 2384 pixels and 50 spectral bands with wavelengths ranging from 380 to 1050 nm sampled at 10 nm intervals. The scene contains 20 urban landcover classes. The false-color composite image, ground truth image, and class names for the Houston 2018 dataset are provided in Figure 8.

**Figure 8.** (**a**) False-Color Composite Image, (**b**) Ground Truth Image, and (**c**) Class names for the Houston-2018 Dataset.

#### *3.2. Experimental Setup*

Our proposed approach was evaluated by comparing its performance with eight state-of-the-art approaches for HSI feature extraction (Algorithm 2, see Section 3.4.5). These include SVM [7], Edge Preserving Filter (EPF) [26], superpixel-based classification via multiple kernels (SCMK) [34], region-based relaxed multiple kernel (R2MK) [35], adjacent superpixel-based multiscale generalized spatial-spectral kernel (ASMGSSK) [36], Multiscale superpixel-based PCA (MsuperPCA) [37], 2D Singular Spectrum Analysis (2D-SSA) [29], and 2D Multiscale Singular Spectrum Analysis (2D-MSSA) [31]. A common way to measure the efficiency of feature extraction is through the accuracy of the classifier scored by the experiments. As a result, the classification setup must be appropriate with the current state-of-the-art. In light of this, SVMs have demonstrated themselves to be robust and efficient in multi-class classification applications.

The LIBSVM toolbox [48] is used to implement SVM as the default classifier for all of the involved methods. A Gaussian RBF kernel is utilized for SVM implementation, and a grid search is applied to tune both key parameters of RBF-SVM; the penalty *c* and the gamma *γ*. The SVM parameters are kept constant across all competitive experiments for a fair comparison. To avoid systematic errors and reduce random discrepancies, all experiments were independently carried out ten times each with different training and testing subsets, with no overlap between each training and the corresponding testing subset. This was intended to ensure good statistical significance for our experiments.

Stratified sampling was used to randomly obtain the training and testing subsets. For training, 3%, 2%, 1%, and 0.2% samples per class were selected for the IP, PU, SAL, and Houston 2018 datasets, respectively. Additionally, four objective quality indices are utilized to evaluate image classification results: namely the OA, the average accuracy (AA), the kappa coefficient, and class-by-class accuracy. All experiments were conducted using MATLAB R2018b software, installed on a personal computer with an Intel core i5-6200 CPU clocked at 2.30 GHz, and 16 GB RAM.

#### *3.3. Parameter Sensitivity Analysis*

Table 1 displays the best parameter settings for the competing algorithms, found by experimentation. For the proposed SP-SSA algorithm, the size of the 2-D embedding window was set to 5 × 5 pixels for the IP and Salinas dataset; whereas, for the PU and Houston 2018 dataset, the window size was set at 3 × 3 pixels. For the IP and SAL datasets, superpixels were set at 100. However, the amount of superpixels in the PU and Houston 2018 datasets were set to 150 and 500, respectively. The effect of window size variation for different number of superpixels on the classification performance for the experimental datasets is provided in Figure 9.

As each superpixel is reconstructed individually, smaller window sizes are preferred since they lead to better image reconstructions. Using a large window may smooth the results too much and result in mixing errors. A 2D-SSA algorithm was presented in [29] for feature extraction in HSI, where various window sizes, such as 5 × 5, 10 × 10, 20 × 20, 40 × 40, and 60 × 60, were examined. The IP and SAL datasets produced the best classification accuracy when the window size was set at 10 × 10. When analyzing the PU and Houston dataset, the window sizes of 5 × 5 showed the best classification results. Since the optimal window size may vary depending on the dataset, ref. [31] adopts a multiscale strategy to improve the generalization ability.

**Figure 9.** Effect of window size variation for different number of superpixels on the classification performance for the (**a**) Indian Pines (**b**) Pavia University, and (**c**) Salinas, and (**d**) Houston 2018 datasets.


**Table 1.** Parameter settings for different algorithms in the Indian Pines, Pavia University, Salinas, and Houston 2018 Datasets.

#### *3.4. Experimental Result and Analysis*

In this section, the four HSI data sets outlined in Section 3.1 are utilized, and several experiments are performed to examine the efficacy of the proposed SP-SSA method. Figure 10 compares the classification results obtained with varying numbers of training samples on four datasets. It can be noted that better classification performance is evident when larger numbers of labeled samples are utilized for training; after passing the percentages used in this work; however, the accuracy level mostly plateaus, and no further significant improvement is observed. Our proposed approach attains the best classification accuracy in almost all cases, regardless of the number of samples, proving its robustness. Classification results from all four data sets are provided in Tables 2–5 and quantitatively support the dominance of the proposed method. Individual classification maps generated by the proposed SP-SSA method and all the compared approaches are displayed in Figures 11–14.

**Figure 10.** Effect of training sample variation on the classification performance for the (**a**) Indian Pines (**b**) Pavia University, (**c**) Salinas, and (**d**) Houston 2018 datasets.

#### 3.4.1. Results from the Indian Pines Dataset

Based on the results shown in Table 2, the proposed method achieves the best values across three metrics, and its accuracy exceeded 89% on almost all classes. In the tables, the best results in each row are highlighted in bold font. When comparing SP-SSA with raw HSI data, the OA improved substantially from 76.42% to 98.15%. In addition, comparisons between SVM and other methods indicated that the incorporation of spatial features can enhance the classification performance compared to considering spectral features alone.

Superpixel-based methods, such as SCMK, R2MK, ASMGSSK, and MsuperPCA techniques, yield higher classification accuracy as compared to non-superpixel based techniques (EPF, 2D-SSA, and 2D-MSSA); by grouping spectrally identical regions, superpixels offer a powerful way to exploiting spatial/contextual information. It can also be noted that methods considering multi-scale windows (ASMGSSK, MsuperPCA, and 2D-MSSA) perform better with respect to fixed-window methods. Due to the different window sizes, unique local spatial features can be exploited, which allows better covering of different sizes of land cover classes and different scales of spatial features. On the downside, the use of multiscale approaches involves heavier processing burdens. As the proposed method reconstructs each superpixel individually, better classification results are obtained.

**Table 2.** Classification results for the Indian Pines Dataset with 3% training for SVM, EPF, SCMK, R2MK, ASMGSSK, MsuperPCA, 2D-SSA, 2D-MSSA, and SP-SSA algorithms.


Figure 11 displays the classification maps produced by various approaches for the Indian Pines dataset. For the SVM approach, the classification map appears very noisy if spatial features are not considered. Through the use of neighborhood spatial information, the EPF and 2D-SSA techniques can suppress spot-wise misclassification to a large extent, but these methods do not preserve the detailed structures of the HSI well enough.

However, by adopting superpixel-based approaches, the generated classification map becomes much smoother, and more accurate estimates are obtained in the detailed region. With the utilization of multi-scale approaches (like ASMGSSK, MSuperPCA, and 2D-MSSA), the amount of misclassification is further reduced. Still, even with multi-scale approaches, landcover boundaries are frequently misplaced. As can be observed from Figure 11, the proposed approach effectively solved the above-mentioned problems due to its considerate utilization of spectral and spatial features.

**Figure 11.** (**a**) Ground Truth Image, Classification Maps of (**b**) SVM (**c**) EPF (**d**) SCMK (**e**) R2MK (**f**) ASMGSSK (**g**) MsuperPCA (**h**) 2D-SSA (**i**) 2D-MSSA (**j**) SP-SSA for Indian Pines dataset.

#### 3.4.2. Results from the Pavia University Dataset

Quantitative results are presented in Table 3. The proposed SP-SSA method still achieved higher classification accuracy and ranked first among all the compared methods, closely followed by the ASMGSSK algorithm. Also, in comparison to EPF, SCMK, R2MK, 2D-SSA, MSuperPCA, and 2D-MSSA techniques, the average improvement of the proposed approach is over 4.41%, 3.64%, 2.09%, 2.37%, 1.3%, and 1.48%, respectively. For comparison, the top results in the tables are boldfaced. In Figure 12, different classification maps are shown, based on various testing methods applied to the PU dataset.

According to Figure 12, the classification map for SVM still continues to remain noisy. Both EPF and 2D-SSA can generate a relatively smooth result; however, some significant regions remain undetected (e.g., the detailed areas). The superpixel-based methods (SCMK, R2MK, ASMGSSK, and MSuperPCA) and SSA-based approach (2D-SSA and 2D-MSSA) offer significantly improved performance, but the proposed 2D-SSA method remains the most promising approach as it outperforms all the compared algorithms.


**Table 3.** Classification results for the Pavia University Dataset with 2% training for SVM, EPF, SCMK, R2MK, ASMGSSK, MsuperPCA, 2D-SSA, 2D-MSSA, and SP-SSA algorithms.

**Figure 12.** (**a**) Ground Truth Image, Classification Maps of (**b**) SVM (**c**) EPF (**d**) SCMK (**e**) R2MK (**f**) ASMGSSK (**g**) MsuperPCA (**h**) 2D-SSA (**i**) 2D-MSSA (**j**) SP-SSA for the Indian Pines dataset.

#### 3.4.3. Results from the Salinas Dataset

The visual classification maps and quantitative results obtained by various classifiers on the Salinas dataset are shown in Figure 13 and Table 4, respectively. In the table, the best results are shown in bold. Based on the visual quality as well as objective metrics, it can be observed that the proposed SP-SSA method outperformed other competing approaches. In addition, compared with the 2D-SSA method that globally reconstructs the image using fixed-size embedded windows, the SP-SSA method considers the local spatial information by reconstructing each superpixel individually, which helps in further reducing the disturbances and improving the class assignment.

**Table 4.** Classification Results for the Salinas Dataset with 1% training for SVM, EPF, SCMK, R2MK, ASMGSSK, MsuperPCA, 2D-SSA, 2D-MSSA, and SP-SSA algorithms


**Figure 13.** (**a**) Ground Truth Image, Classification Maps of (**b**) SVM (**c**) EPF (**d**) SCMK (**e**) R2MK (**f**) ASMGSSK (**g**) MsuperPCA (**h**) 2D-SSA (**i**) 2D-MSSA (**j**) SP-SSA for Salinas dataset.

#### 3.4.4. Results from the Houston 2018 Dataset

The quantitative results for the Houston 2018 dataset with 0.2% training samples from each class are presented in Table 5. The corresponding classification map is shown in Figure 14. The best results from the tables are displayed in bold font for comparison. As observed from Table 5, the proposed methods are robust and achieve good classification results even for challenging scenes. The proposed approach improves accuracy from 68.19% to 83.57% for the SVM method. In this case also, the superpixel-based approaches (SCMK, R2MK) display superior performance as compared to non-superpixel based methods (EPF, 2DSSA). Here also, multi-scale window approaches (ASMGSSK, MsuperPCA, and 2D-MSSA) outperform fixed-window based methods as different scales of spatial features are incorporated into the analysis. Figure 14 also highlights the superiority of the proposed method. The salt and pepper noise is reduced by a greater extent, and a smoother classification map is produced with the proposed method.




**Table 5.** *Cont.*

**Figure 14.** (**a**) Ground Truth Image, Classification Maps of (**b**) SVM (**c**) EPF (**d**) SCMK (**e**) R2MK (**f**) ASMGSSK (**g**) MsuperPCA (**h**) 2D-SSA (**i**) 2D-MSSA (**j**) SP-SSA for Houston 2018 dataset.

#### 3.4.5. Statistical Evaluation

The effectiveness of the proposed method was statistically evaluated using McNemar's test. The classification results for all the test cases were compared using this test. The McNemar's test is defined as in Equation (12), where it is assumed that two generic algorithms, named Algorithm 1 and Algorithm 2 are compared.

$$Z = \frac{f\_{12} - f\_{21}}{\sqrt{f\_{12} + f\_{21}}} \tag{12}$$

In the equation above, *f*<sup>12</sup> indicates the number of samples correctly classified by Algorithm 1 and incorrectly classified by Algorithm 2, and *f*<sup>12</sup> indicates the number of samples for the opposite case. The performance of Algorithm 1 is better than Algorithm 2 if *Z* > 0. The differences between Algorithm 1 and Algorithm 2 are statistically significant if |*Z*| > 1.96. In our case, Algorithm 1 is the algorithm proposed in our manuscript, and Algorithm 2 is —sequentially— each one from the list of standard algorithms: SVM, EPF, SCMK, R2MK, ASMGSSK, MsuperPCA, 2DSSA, 2DMSSA.

McNemar's test between the proposed SP-SSA algorithm and the algorithms listed above for the Indian Pines, Pavia University, Salinas, and Houston 2018 datasets are provided in Table 6. The test result clearly reveals that the classification results for the

proposed method were significantly better—in a McNemar's statistical sense—compared with other approaches.


**Table 6.** Statistics of the McNemar Test for the Indian Pines, Pavia University, Salinas, and Houston 2018 datasets.

#### *3.5. Advantage of Proposed Method over 2D-SSA*

#### 3.5.1. Applying SP-SSA on General Images

In the proposed approach, 2D-SSA is applied on each and every superpixel segmented region. Hence, it can be considered as a local 2D-SSA approach that can extract accurate spatial information on each single object. In the case of global 2D-SSA, features are oversmoothed, and features are not prominent for specific classes. In local 2D-SSA instead, object-specific texture information can be highlighted. In Figure 15, the popular cameraman image and an artificial test image are used to demonstrate the effectiveness of the proposed approach over the 2D-SSA approach.

When the cameraman image is reconstructed using the 2D-SSA method, the Mean Square Error (MSE) comes out to 115.8865; however, when the same image is reconstructed using the proposed SP-SSA approach, the MSE reduces to 93.0468. A similar behavior is also observed with the test image. With the proposed SP-SSA method, the MSE reduces to 237.1038 from 287.5323. This signifies that the proposed method can reconstruct an image with minimum error and can effectively integrate local information during the reconstruction process.

**Figure 15.** (**a**) Cameraman image (**b**) 2D-SSA Reconstructed image [MSE = 115.8865] (**c**) SP-SSA reconstructed image [MSE = 93.0468] (**d**) Test image (**e**) 2D-SSA Reconstructed image [MSE = 287.5323] (**f**) SP-SSA reconstructed image [MSE = 237.1038].

#### 3.5.2. Applying SP-SSA on HSI

The HSI is composed of a stack of 2D images carrying valuable information about each spectral band. To demonstrate the effectiveness of the proposed method, a randomly selected spectral band at 667 nm was considered for our analysis. Figure 16b,c contains the scene as reconstructed by 2*D*-*SSA* and *SP*-*SSA*, respectively. Since the HSI was acquired over a large area, it includes multiple objects with different textural information. This is a typical case where object-specific reconstruction works better than direct reconstruction. Textural information can be highlighted accurately by using local reconstruction as opposed to global reconstruction. The error in SP-SSA-based reconstruction is indeed lower as compared to 2D-SSA-based reconstruction. The same conclusion can also be drawn from Figure 16.

In the case of 2D-SSA-based reconstruction, the Mean Square Error (MSE) is 612.4349, while, in the case of SP-SSA-based reconstruction, the MSE is 504.5685. Figure 16d,e contains the difference image for 2D-SSA-based reconstruction and SP-SSA-based methods. It can be clearly observed that edge information is preserved with the proposed method. The SP-SSA-based reconstruction is applied to all spectral bands and generates a modified hypercube with preserved local structure information and minimum noise level. These latter features generally lead to better classification performance.

**Figure 16.** (**a**) Original scene at band 667 nm (**b**) Reconstructed scene by 2D-SSA (**c**) Reconstructed scene by SP-SSA (**d**) Difference image for the 2D-SSA reconstructed scene (**e**) Difference image for the SP-SSA reconstructed scene.

#### **4. Conclusions and Future Scope**

Feature extraction is one of the most crucial steps in HSI classification. It is essential to capture comprehensive spatial and spectral information for accurate feature extraction. For image reconstruction, the conventional 2D-SSA algorithm usually extracts spatial features directly by applying the embedding window to the entire image. However, HSI scenes frequently encompass a broader area and contain several items. As a result, spatial information pertaining to local objects must be recovered. To solve this problem, in the proposed method, a superpixel-based SSA technique was presented, which can capture the object specific spatio-spectral information accurately.

In this work, the original HSI was first divided into various semantic sub-regions by the superpixel segmentation algorithm. Next, each segment was reconstructed individually by applying 2D-SSA. The generated reconstructed HSI was then classified using the SVM classifier, and the final classification map was produced. Local characteristics may be collected effectively in the suggested method since 2D-SSA is applied at the superpixel level. However, two parameters must be adjusted: the amount of superpixels and the embedding window size. Future developments will aim at finding the optimal criteria to determine the parameters of the procedure and to investigate relationships between the characteristics of the HSI and quality of the results.

**Author Contributions:** Conceptualization, P.K.B. and R.P.; methodology, S.S.; software, R.P. and S.S.; validation, S.S., R.P. and P.K.B.; formal analysis, S.S.; investigation, S.S. and R.P.; resources, F.D. and P.K.B.; data curation, S.S.; writing—original draft preparation, S.S.; writing—review and editing, F.D.; visualization, S.S. and F.D.; supervision, F.D.; project administration, F.D. and P.K.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available datasets have been utilized for conducting the experiments. Details are provided in the text of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Review* **Spatial Analysis for Landscape Changes: A Bibliometric Review**

**Maria Danese \* and Dario Gioia**

Institute of Cultural Heritage-National Council of Research, C.da S. Loja, Tito Scalo, 85100 Potenza, Italy; dario.gioia@cnr.it

**\*** Correspondence: maria.danese@cnr.it; Tel.: +39-0971427324

**Abstract:** The main aim of this study is to analyze from a bibliometric point of view the research trend in spatial analysis for landscape changes using the records published in the Web of Science database in the last twenty years. Several parameters such as documents published per year, sources of documents, number of citations as well as VOSviewer software and GIS are used for the analysis of different metrics such as the number of citations, co-authorship network, and keyword occurrences. Analysis of the number of papers, their keywords, and authorships countries shows the research trend in the specific topics of the spatial analysis for landscape changes and consequently can constitute a benchmark for researchers who approach this research topic.

**Keywords:** spatial analysis; landscape changes; bibliometric mapping; Web of Science

#### **1. Introduction**

Spatial analysis saw a big improvement in the last few decades thanks to the parallel growth of Information Computer Technologies (ICT), both from the hardware and from the software development point of view. Since then, many research fields have been taking advantage of the spatial analysis discipline. Between these, there are all the studies directly concerned with human- and natural-induced landscape changes. Bibliometric analysis of "hot" research topics is a growing tool for the investigation of emerging disciplines, cooperations, publication impact, and new research trends [1]. In recent years, such an approach has been adopted by several authors to quantitatively delineate the global trend of different research topics or wider disciplines [2–6], highlighting the usefulness of the method to study the past and future direction of research patterns. Moreover, the recent availability of wider databases and innovative software of bibliometric analysis provides new tools for the deeper visual and statistical analysis of temporal and geographic global distribution in research trends [1]. In this work, we carry out a bibliometric analysis of the emerging research topic of the spatial analysis of landscape changes with the aim to investigate the topic's research pattern in the last two decades and to guide researchers to understand future trends.

Landscape changes are very important for the past, present, and future of the Earth and consequently for human life; consequently, as it is shown in this paper, the study of landscape changes through the help of spatial analysis has grown in recent decades and become a hot topic. For this reason, there are many articles belonging to the Ecology, Territorial Planning and Earth Science sectors. Just some examples of the themes afforded in these areas are papers that use spatial analysis for understanding how landscape changes impact species distribution over space and their interaction with environments [7,8], the quantification of natural resources, such as water [9], forests [10–12], and, more in general, of ecosystem services [13,14], and multi-temporal analysis of satellite images and maps as a tool to reconstruct land-use changes and rates of geomorphological processes [15–17]. One of the most effective approaches to investigate natural landscape changes is the geomorphological one, which benefits from recent advances in the development of software

**Citation:** Danese, M.; Gioia, D. Spatial Analysis for Landscape Changes: A Bibliometric Review. *Appl. Sci.* **2021**, *11*, 10078. https:// doi.org/10.3390/app112110078

Academic Editor: Nadia Giuffrida

Received: 24 September 2021 Accepted: 25 October 2021 Published: 27 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and algorithms of digital elevation model (DEM) comparison [18,19], image change detection [20,21], and landscape evolution models [22–25]. Moreover, sustainability is a big topic connected with this research line and there are many articles concerning it [26–28], together with analysis connected to climate change [29–31] and pollution quantification [32–34]. However, there are other disciplines such as Archaeology or Historical studies [35–37], or fire analysis [38–40], that are strictly connected to the wider topic of landscape changes.

From what has been said, it is now clear that to consider "spatial analysis for landscape changes" (from here s.a.l.c.) implies covering a very broad area because there are many topics and, at the same time, there are many types of spatial analysis that could be used to investigate these changes, from more traditional methods [41] to more intelligent data analysis [42].

Consequently, the main aim of this paper is to clarify topics, trends, and methods that are connected to the research line of s.a.l.c. through a comprehensive and accurate bibliometric analysis during the last twenty years.

#### **2. Methodology**

#### *2.1. Bibliometric Search Engine, Tools, and Software Used for S.A.L.C. Analysis*

To retrieve the scientific literature for this paper, we used the Web of Science database (WoS). Among the most popular bibliographic databases of research papers such as Scopus, WoS, and Google Scholar, we chose WoS because it represents the best compromise among database completeness and cataloguing of higher rank journals [43].

This research was conducted by using the combination of s.a.l.c. keywords, written between the quotation marks, in order to include at the same time all the four keywords and by looking for them inside the title, abstract, and keywords of the scientific paper. Then, the analysis excluded from the research all the WoS categories including the life science disciplines, such as biotechnology, applied microbiology, immunology, oncology, etc.

Regarding the time of publication, two temporal ranges were considered.

The first range considers years from 2001 to 2020, in order to study the general characteristics of scientific literature in the highlighted topic.

The second range considers the same period, from 2001 to 2020, but is divided into four five-year intervals, in order to investigate the long-term literature trends.

To achieve this goal, the following tools were used: first, some tools present inside WoS were used to extract some statistics and graphs; second, an open-source Geographic Information System, QGIS (open source software, downloadable at the site https://www. qgis.org/it/site/, accessed on 10 October 2021), was used to analyze the geographic distribution of selected works; third, for more specific diagrams, the VOSviewer software was used. It is a free software [44] based on text-mining and it is useful to construct bibliometric maps.

The bibliometric mapping discipline, also called science mapping, recently appeared in order to look for relationships between documents, keywords, or authors [45,46], and with the main aim to investigate the structure and the dynamics of a topic [1,47].

In the literature, other software programs exist for bibliometric mapping, such as SCImat [48] or CiteSpace [49]. Some software reviews are present in [50,51]. We chose VOSviewer for its simplicity, flexibility, and clearness of the results.

#### *2.2. VOSviewver Diagrams*

Concerning the use of VOSviewer, it offers many useful tools for bibliographic analysis. In this paper, the following diagrams were used:

• The keyword co-occurrence map. It is a distance graph showing the connection between the keywords included in the selected bibliography. If the terms co-occur inside the same phrases, a higher relevance score is assigned to them. Consequently, terms that are linked and near each other in the map are more related. With the co-occurrence map, it is then possible to analyze the main keywords that characterize the state of the art in a domain field. With this type of analysis, VOSviewer returns the

graph showing links between key terms, which are also divided into clusters; these clusters are in turn built up on the basis of the co-occurrences of terms inside the paper titles. Moreover, the number of occurrences and the total link strength are associated with each term.


#### **3. Results**

#### *3.1. General Quantitative Results*

The first result, obtained from the research in Web of Science, is the quantification of the citations and publications existing in the research field s.a.l.c. A total of 8409 records were found, with a minimum of 86 in 2001 and a maximum of 862 in 2020. Figure 1 shows a consistent increase both in the publication and in the citation number, except for two light drops between years 2002–2003 and 2016–2017. There was also a little decrease in 2020. However, we have to remember that 2020 was an exceptional year because of the COVID-19 pandemic, and this could influence the number and the topic of the worldwide research papers.

**Figure 1.** Diagram representing the number of citations and publications reported by Web of Science from 2000 to 2020.

#### *3.2. Result Heterogeneity*

In the results obtained by Web of Science, the first aspect that emerges, as expected and explained in the introduction, is the heterogeneity of categories. This aspect is clearly highlighted in Figure 2, where the composition of the first ten categories is visualized in a tree map. Visual inspection of the map highlights that the main categories working in these sectors are environmental disciplines such as Environmental Sciences (in 1st place with 2349 documents, corresponding to 27.6% of the total data), Ecology (2nd place, 2337, 27.4%), Multidisciplinary Geosciences (3rd place, 1323, 15.5%), and Physical Geography (4th place, 1236, 14.5%). Moreover, it is possible to discover a first trend, not directly asked with the query conducted in the Web of Science database. In fact, categories such as Remote Sensing (873, 10.2%) and the Imaging Science Photographic Technology categories (528, 6.2%) are, respectively, classified in 5th and 9th place, suggesting the first delineation of the main techniques and approaches related to the s.a.l.c.

**Figure 2.** Tree map of the first twenty Web of Science categories where the works about spatial analysis for landscape archaeology are published.

#### *3.3. Spatial Distribution of S.A.L.C.*

S.a.l.c. literature authorship is spatially distributed all around the world, even if a preliminary analysis of the results indicates a relevant prevalence of European-related papers coming from 129 different countries. Figure 3 shows the frequency of documents found for each country and classified in deciles. The first five countries interested in the research area are the United States of America (with 2765 documents, 32.9%), China (1205, 14.3%), Germany (888, 10.0%), the United Kingdom (836, 10.0%), and Canada (628, 7.5%). It is worth noting that the countries that seem not involved in the research topics and present zero articles or very few articles (1 or 2) are prevalently located in the African continent.

It is also interesting to look at the co-authorship map between countries (Figure 4), and the four main geographical clusters found that are those with the following rank: (1st) the USA, (2nd) China, (3rd) Germany, and (4th) the Netherlands.

#### *3.4. Co-Occurrence Map of Keywords of the Whole Period*

Concerning the co-occurrence map calculated in the period between 2016 and 2020, four main clusters were found. The prevailing keywords of these clusters are the following (the number of the cluster is just nominal, not ordinal):


instead of the "method" cluster (the Cluster 3), even if, of course, the link with it remains through keywords GIS, land-use change, remote sensing, and, most of all, urbanization, which is strictly linked with topics such as the protection or conservation of biodiversity.


**Figure 3.** Spatial distribution of the countries of article authorship (**a**) in the World and (**b**) in Europe. World background taken from @ naturalearthdata.com. The graph was extracted using QGIS software.

**Figure 4.** Co-authorship cluster map. Graph obtained by VOSviewer software.

(**a**) **Figure 5.** *Cont*.

**Figure 5.** (**a**) Co-occurrence map (period: 2001–2020); (**b**) overlay visualization between key terms and their citation year. Graphs made with VOSviewer software.

By "main word", we meant words with the major frequencies clustered by the software. Together with the co-occurrence map, we extracted a map where the overlay between terms and the citation year is shown (Figure 5b): it shows the most cited keywords over the years. In the legend, in particular, the years between 2012 and 2016 are highlighted because this is the temporal range with the largest number of citations. In particular, the most recently cited keywords since 2016 are climate change and ecosystem service, while "older" ones (and actually, they are still consolidated in the literature) are GIS and remote sensing.

#### *3.5. Density Map for Interval of 5 Years*

In order to investigate the trend in s.a.l.c. research topics on a shorter-term period, we extracted the density map of keywords for four different intervals of 5 years. The interpretation of the maps (Figure 6) could become clearer with the help of data and histograms reported in Figure 7, where occurrence frequencies of the keywords are reported and schematized. Visual inspection of the maps suggests that the two common keywords in the top five list are "landscape" and "pattern". These two words are clearly connected to the main topic of this research. Instead, we can see, for the first two temporal blocks (2001–2005 and 2006–2010), that the other keywords are "model" and "GIS", which are related to the methodological aspects of the research. These two terms, in particular, have a great increase in the second time block (i.e., from 2006 to 2010). Here, the keyword "climate change" also appears. The latter becomes more significant in the third block (2011–2015) and also increases in the 2016–2020 interval. The keyword "climate change" is further reinforced by the last keyword of the top five list, respectively, in 2011–2015, with the terms "biodiversity" and "conservation", and in 2016–2020, with the terms "impact" and, again, "conservation".

Finally, Figure 7 shows the number of occurrences and the related histograms of the keywords already cited in Section 3.4, which considers the most important terms for each

cluster found in the total period between 2001 and 2020. This is performed to observe the trend of main keywords (in bold, Figure 7) and some of the secondary keywords.

**Figure 6.** Density map. Obtained with VOSviewer.


**Figure 7.** Occurrences and related histograms for the main keywords observed in the 2001–2020 period in order to observe their trend in the four periods.

#### *3.6. Types of Spatial Analysis*

Although spatial analysis is a relatively broad term adopted for different and wide research disciplines, an analysis of keywords can be useful to extract information about the main methods and approaches used in s.a.l.c research. To this aim, data were cleaned manually and only keywords concerning methods were left in the analysis. Moreover, keywords with a broader meaning such as "spatial analysis", "remote sensing", or "geospatial analysis" were excluded.

A synoptic scheme of such an analysis is reported in Figure 8, where one can observe the overlay between the co-occurrences map and its overlay with the publication year. Table 1 shows the first 50 keywords, in order of importance for their occurrences: the two top keywords represent, respectively, the traditional method used to analyze a landscape and one of the more innovative methods. In fact, the keyword "classification" (359 occ.) can be correlated to the results of the landscape analysis, whereas the second one (i.e., "simulation", 205 occ.) has a strong connection with the approaches adopted for the validation of the results or landscape modeling. Other relevant keywords are the tools used by spatial analysis, such as "indicators" (170 occ.) and "metrics" (156 occ.), or more innovative and specific methods such as "cellular automaton" (113 occ.), "change detection" (104 occ.), "spatial autocorrelation" (103 occ.), "gradient analysis" (84 occ.), "species distribution model" (82 occ.), and "regression" (80 occ.). Among them, "cellular automation" is maybe the most innovative, since has been recently introduced and represents a spatially distributed evolution of artificial intelligence algorithms. Finally, in 10th place, we find the term "(DEM, 72 occ.)", which represents the basic elements of most of the landscape analysis approaches. Additional information about research trends in the s.a.l.c. topic can be inferred from the temporal variation of the keyword occurrences: for example, the prevalence of keywords such as DEM, statistical analysis, data analysis, or regression in the 2012–2013 period (blue tones in Figure 8) seems to suggest research approaches mainly based on visual inspection and basic statistical analysis, while the appearance of peculiar terms in 2016 such as random forest, Markov chain, machine learning, or cellular automata clearly indicates a transition toward the automatic or semi-automatic classification of landscape changes.

**Figure 8.** Overlay map for the types of spatial analysis more occurrent in s.a.l.c. research topic. Made with VOSviewer.


**Table 1.** First 100 keywords with the highest number of occurrences.

<sup>1</sup> Occurrences.

#### **4. Discussion and Concluding Remarks**

In this paper, a detailed bibliometric investigation of the research trend in spatial analysis for landscape changes using the records published in the Web of Science database in the last twenty years was carried out. Such an approach has been widely used to drive scientists that need to understand the hot topic of a research field or territory (see, for example, [2–4,52–54]).

Our analysis was conducted with the help of three different software packages: the analysis utility offered by WoS, a GIS software (QGIS), and a specific software for bibliometric mapping, VOSviewer. The results highlight that the topic has received increasing attention in the last two decades. As a matter of fact, we observed a constant and exponential increase in the number of papers and citations since 2000. Such an increase can be partly ascribed to: (i) the growing availability of high-resolution DEMs and remote sensing images; (ii) automatic tools or algorithms of landscape classification and the analysis of land-use changes.

Our results also suggest that the research topic is multidisciplinary, ranging from different disciplines such as Environmental Sciences, Earth Sciences, and Ecology, and it is mainly conducted by Chinese, European, and North American research groups. Moreover, on the basis of the statistical analysis of keyword occurrences, it is possible to reconstruct the following main research patterns:


To achieve a deeper understanding of the research trend, each map in a wide research field such as s.a.l.c. would request a widening of the analysis and the reading of results in the different clusters and thematic areas identified. This would provide enough material for many other future studies in order to identify an internal state of the art and trends of those subsectors.

**Author Contributions:** Conceptualization, M.D. and D.G.; methodology, M.D. and D.G.; software, M.D and D.G.; validation, M.D. and D.G.; formal analysis, M.D. and D.G.; investigation, M.D. and D.G.; resources, M.D. and D.G.; data curation, M.D. and D.G.; writing—original draft preparation, M.D. and D.G.; writing—review and editing, M.D. and D.G.; visualization, M.D. and D.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the data used in this paper can be found at the Web of Science search engine (https://www.webofscience.com/wos/woscc/basic-search, accessed on 10 October 2021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com