Exploring Interpretability in Deep Learning-Based Forest Classification - Separating Similar Forest Classes using Sentinel-2 Time Series

# Exploring Interpretability in Deep Learning-Based Forest Classification
Source: [https://arxiv.org/html/2610.19270](https://arxiv.org/html/2610.19270)
firstpage,color=gray\!60,angle=90,scale=0\.32, xpos=\-4\.05in,ypos=0\][Publication doi](https://doi.org/)firstpage,color=gray\!60,angle=90,scale=0\.32, xpos=3\.9in,ypos=0\][Preprint doi](https://doi.org/)firstpage,color=gray\!90,angle=0,scale=0\.28, xpos=0in,ypos=\-5in\]\*correspondence:maximilian\.kirsch@h\-ka\.de

## Exploring Interpretability in Deep Learning\-Based Forest Classification \- Separating Similar Forest Classes using Sentinel\-2 Time Series

Maximilian Kirsch[ID](https://orcid.org/0009-0001-1754-7809)††thanks:maximilian\.kirsch@h\-ka\.deAffiliation:Intelligent Systems Research Group \(ISRG\), Karlsruhe University of Applied Sciences, GermanyChristine Preisach[ID](https://orcid.org/0009-0009-1385-0585)Affiliation:Intelligent Systems Research Group \(ISRG\), Karlsruhe University of Applied Sciences, Germany

###### Abstract

The ecological insurance theory posits that diverse forests are more resilient and better able to adapt to changing environmental conditions\. Increasing drought and heat stress in Central Europe promotes mixed forests for maintaining ecosystem services\. However, a diverse species distribution complicates remote sensing\-based species mapping, as distinguishing pure from mixed species compositions at the pixel level remains challenging\. We evaluated the potential of recurrent neural networks \(RNNs\) for tree\-species classification from Sentinel\-2 time series and investigated the spectral\-temporal features underlying class discrimination, with particular emphasis on separating pure from mixed Norway spruce pixels\. We compared RNN classification models with Random Forest \(RF\) and Extreme Gradient Boosting \(XGBoost\) and used Integrated Gradients \(IG\) to quantify feature attributions across spectral bands and time\. The best\-performing RNN achieved a Macro\-F1 score of 0\.851 on the test data, outperforming XGBoost by 0\.012 and Random Forest by 0\.103\. Beyond classification performance, the attribution analysis revealed distinct spectral\-temporal patterns associated with the separation of pure and mixed Norway spruce pixels\. Attribution values increased particularly after approximately day of year \(DOY\) 125, while B11 showed the strongest contribution during periods of increasing spectral divergence between the two classes\. These findings show that recurrent models can provide competitive classification performance while also revealing the phenological periods and spectral regions that are most relevant for distinguishing closely related forest compositions\. Combining time\-series classification with explainability therefore provides additional insight into the spectral\-temporal mechanisms underlying tree\-species discrimination and may support more detailed mapping of mixed forest structures\.

Preprint, compiled October 13, 2026

*K*eywordsTree Species Classification⋅\\cdotSITS⋅\\cdotDeep\-Learning⋅\\cdotExplainability⋅\\cdotXAI

## 1Introduction

Understanding the spatial distribution of tree species provides a basis for decision\-making in forest management, biodiversity conservation, and climate protection\[[1](https://arxiv.org/html/2610.19270#bib.bib21)\]\. Repeated monitoring of forest composition can further reveal changes caused by natural disturbances, harvesting, and regrowth, thereby supporting forest inventories and long\-term forest management\[[2](https://arxiv.org/html/2610.19270#bib.bib13),[3](https://arxiv.org/html/2610.19270#bib.bib22)\]\. Under increasingly frequent heat, drought, and disturbance events\[[4](https://arxiv.org/html/2610.19270#bib.bib3)\], forests in Central Europe need to be more resilient and adapted to change\[[5](https://arxiv.org/html/2610.19270#bib.bib10)\]\. The ecological insurance theory posits that greater species diversity can stabilize ecosystem functioning to environmental fluctuations\[[6](https://arxiv.org/html/2610.19270#bib.bib5)\]\. Hence, theNew EU Forest Strategy for 2030of the European Commission explicitly promotes the creation and maintenance of functionally diverse, mixed\-species forests\[[7](https://arxiv.org/html/2610.19270#bib.bib41)\]\. Consequently, mapping approaches must also reliably represent mixed species compositions\. Satellite\-based forest classification has been applied across a wide range of spatial scales, from site\-level studies to global forest mapping\[[8](https://arxiv.org/html/2610.19270#bib.bib12),[9](https://arxiv.org/html/2610.19270#bib.bib11)\]\. Sentinel\-2 is particularly well suited for this purpose because it provides freely available multispectral imagery at high revisit frequency and spatial resolutions of up to 10 m\. Early applications largely relied on single\-date imagery\. For example,[Immitzer et al\. \[10\]](https://arxiv.org/html/2610.19270#bib.bib19)evaluated Sentinel\-2 data for crop and tree species classification in Central Europe using a supervisedRandom Forest \(RF\)classifier and achieved an overall accuracy of 65% for tree species classification\. Their results further highlighted the importance of the red\-edge and short\-wave infrared \(SWIR\) bands\. However, single\-date observations capture only one stage of the seasonal development and may therefore provide limited discriminatory information for species with similar spectral characteristics at a given phenological stage\. With the launch of Sentinel\-2B in 2017, the increased observation frequency facilitated the use ofSatellite Image Time Series \(SITS\)for vegetation and forest classification\[[11](https://arxiv.org/html/2610.19270#bib.bib15)\]\. Exploiting this temporal information can improve species discrimination by capturing phenological differences in seasonal spectral trajectories\.[Grabska et al\. \[12\]](https://arxiv.org/html/2610.19270#bib.bib20), for example, compared single\-date and multi\-temporal Sentinel\-2 imagery for the classification of nine tree species in mixed forests and reported an increase in overall mapping accuracy of 5–10% when using time\-series data, with combinations of spring and autumn observations being particularly informative\. In tree species classification, Sentinel\-2SITSare commonly analyzed using machine\-learning approaches, with tree\-based ensemble methods such asRFandExtreme Gradient Boosting \(XGBoost\)being frequently applied \(e\.g\.,\[[12](https://arxiv.org/html/2610.19270#bib.bib20),[13](https://arxiv.org/html/2610.19270#bib.bib17),[14](https://arxiv.org/html/2610.19270#bib.bib33),[15](https://arxiv.org/html/2610.19270#bib.bib31),[16](https://arxiv.org/html/2610.19270#bib.bib32)\]\)\. These methods are widely used because of their strong predictive performance, scalability to large study areas, and comparatively high interpretability\[[17](https://arxiv.org/html/2610.19270#bib.bib16)\]\. However, SITS contain sequential spectral information whose temporal dependencies can also be modeled directly using deep learning approaches\. The continued operation of Sentinel\-2 has progressively expanded the available image archive, providing increasingly long time series that also facilitate the application of such sequence\-based deep\-learning approaches \(e\.g\.,\[[18](https://arxiv.org/html/2610.19270#bib.bib14),[17](https://arxiv.org/html/2610.19270#bib.bib16)\]\)\.

Despite their increasing relevance for ecology and management, distinguishing mixed forests from pure stands using remote sensing data is challenging\[[19](https://arxiv.org/html/2610.19270#bib.bib39)\]\. In medium\-resolution satellite images, such as Sentinel\-2 \(10m\), pixels are likely to contain several individual trees, and the proportion of each species within the pixels can vary substantially, which is reflected in classification performance\.[Blickensdörfer et al\. \[13\]](https://arxiv.org/html/2610.19270#bib.bib17), for example, mapped the main tree species in Germany using combined Sentinel\-1 and Sentinel\-2 time series and integrated exogenous features such as temperature and topographic features\. While they reported high F1 scores for the five dominant species in pure stands, ranging from 72% to 97%, the inclusion of mixed stands reduced accuracy by 4\-14 percentage points\. At the same time, most forestry inventory, planning, and management tools have been developed for monocultures\[[5](https://arxiv.org/html/2610.19270#bib.bib10)\]and tree species classification studies prominently focus on pure stands \(e\.g\.,\[[14](https://arxiv.org/html/2610.19270#bib.bib33),[20](https://arxiv.org/html/2610.19270#bib.bib36)\]\)\. This is a limitation, especially in Germany, as mixed forests account for 79% of Germany’s forest area in 2022\[[21](https://arxiv.org/html/2610.19270#bib.bib4)\]\. Disregarding mixed stands in forest mapping products can therefore substantially misrepresent forest composition\. Given the difficulty of distinguishing mixed from pure forest classes, classification performance alone provides limited insight into which spectral and temporal characteristics enable their separation\. Research on explainability in remote sensing\-based tree species classification has predominantly focused on tree\-based methods, partly for their comparatively high interpretability \(e\.g\.,\[[14](https://arxiv.org/html/2610.19270#bib.bib33),[22](https://arxiv.org/html/2610.19270#bib.bib18),[15](https://arxiv.org/html/2610.19270#bib.bib31)\];\[[17](https://arxiv.org/html/2610.19270#bib.bib16)\]\)\.[Immitzer et al\. \[22\]](https://arxiv.org/html/2610.19270#bib.bib18)investigated the spectral variables and acquisition dates most relevant for classifying individual tree species usingRFmodels trained on combinations of Sentinel\-2 bands and spectral indices across 18 scenes\. Their results highlighted the red and short\-wave infrared bands, while observations acquired between April and August contributed most strongly to classification performance\. Similarly,[Wegler et al\. \[14\]](https://arxiv.org/html/2610.19270#bib.bib33)assessed feature contributions inRFandXGBoostmodels using SHAP values and model\-based feature importance and identified the short\-wave infrared and red\-edge bands among the most influential features\. While[Malkanthi et al\. \[23\]](https://arxiv.org/html/2610.19270#bib.bib30)also addressed explainability in deep learning\-based tree species classification, their analysis focused on exogenous variables rather than on the spectral\-temporal information contained in the SITS\. Consequently, although previous studies have provided valuable insights into the spectral and temporal features relevant for tree species classification, the spectral\-temporal decision processes of deep learning models remain insufficiently understood\. This limitation is particularly relevant for closely related classes, such as pure and mixed Norway spruce, for which understanding the wavelengths and phenological periods driving class separation may provide insights beyond overall classification performance\.

To address this gap inSITS\-based tree species classification, we propose an explainability workflow for multiclass classification that identifies the spectral\-temporal features most relevant for class discrimination while remaining applicable to deep learning models\. We particularly investigate the separation of pure and mixed Norway spruce pixels and additionally examine deciduous tree species with similar phenological patterns\. The study addresses the following research questions:

RQ1:How do recurrent deep learning models compare with established tree\-based ensemble models in classifying tree species and distinguishing pure from mixed forest compositions using one\-year Sentinel\-2 time series?

RQ2:Can attribution techniques applied to deep learning models provide insight into the spectral bands and phenological periods relevant for class discrimination, particularly between closely related forest classes?

## 2Materials and Methods

### 2\.1Dataset and Study Area

We used the publicly available Sentinel\-2 tree species time\-series dataset for Germany, published in 2025 by[Freudenberg et al\. \[20\]](https://arxiv.org/html/2610.19270#bib.bib36)\. The dataset combines bottom\-of\-atmosphere reflectance time series from analysis\-ready Sentinel\-2 imagery with tree species labels derived from the German National Forest Inventory\. It covers individual tree observations across Germany from 2015 to October 2022 and includes 48 tree species and three species groups\. We selected the four dominant tree species or species groups in Germany: Norway spruce \(Picea abies\), Scots pine \(Pinus sylvestris\), European beech \(Fagus sylvatica\), and oak, with the latter comprising sessile oak \(Quercus petraea\) and English oak \(Quercus robur\)\[[24](https://arxiv.org/html/2610.19270#bib.bib2)\]\. Furthermore, we split the Norway spruce class into healthy and disturbed Norway spruce, as well as mixed pixels that comprise Norway spruce and other tree species \(mainly beech\) in mixed stands\. Additional classes were deliberately excluded, as the remaining tree species in the dataset by[Freudenberg et al\. \[20\]](https://arxiv.org/html/2610.19270#bib.bib36)were underrepresented for the study’s specific use case of classification and separability analysis\. An overview of the data basis and the spatial distribution of the dataset locations is shown in Figure[1](https://arxiv.org/html/2610.19270#S2.F1)\.

The selected dataset comprised pure and mixed Norway spruce time series, as our objective was to determine whether the models could differentiate between classes that are very closely related\. Most tree\-species maps in Germany focus on pure stands because mixed pixels introduce substantial classification uncertainty\. Their spectral signatures tend to resemble those of the dominant tree species within the pixel, which can lead to mixed pixels being misclassified as pure stands \(see Figure[2](https://arxiv.org/html/2610.19270#S2.F2)\)\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/figs/map_species3_tr_te2.png)Figure 1:Data distribution: spatial distribution of sample locations used in this study, grouped by species\.We subsequently filtered the selected data and applied the preprocessing steps described in Section[3\.1](https://arxiv.org/html/2610.19270#S3.SS1)\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/figs/mixed_tree_classes_copy4.png)Figure 2:Mixed Sentinel\-2 pixels\.Illustration of Sentinel\-2 pixels containing pure Norway spruce stands and mixed stands\.Ashows an orthophoto of a forest structure including Norway spruce and varying deciduous trees\.Bshows a subset of the same area including the outlines of deciduous trees within a 3×\\times3 Sentinel\-2 pixel grid\.Cshows the corresponding Sentinel\-2 pixels averaged over one year \(2021\)\.
### 2\.2Models

In this study, we selected two ensemble methods,XGBoostandRF, as well as four deep learning approaches\.XGBoostandRFwere selected because they have proven to be very successful in the field of SITS forest classification\[[14](https://arxiv.org/html/2610.19270#bib.bib33),[15](https://arxiv.org/html/2610.19270#bib.bib31),[13](https://arxiv.org/html/2610.19270#bib.bib17)\]\. For the deep learning models, we developed an architecture that utilizesRecurrent Neural Networks \(RNNs\), which are particularly well suited to sequential data because they can model temporal dependencies within multivariate time series\[[25](https://arxiv.org/html/2610.19270#bib.bib35)\]\. Their potential forSITSdata has also been explored and validated \(e\.g\.,\[[26](https://arxiv.org/html/2610.19270#bib.bib34),[27](https://arxiv.org/html/2610.19270#bib.bib38)\]\)\. We selectedLong Short\-Term Memory \(LSTM\)\[[28](https://arxiv.org/html/2610.19270#bib.bib40)\]andGated Recurrent Unit \(GRU\)\[[29](https://arxiv.org/html/2610.19270#bib.bib8)\]to extract spectral\-temporal information from Sentinel\-2 time series and embed it into a latent feature space\. Their schematic functionality is illustrated in Figure[3](https://arxiv.org/html/2610.19270#S2.F3)\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/figs/LSTMuGRU_biggerVertical.png)Figure 3:LSTM and GRU architectures\.At time steptt,xtx\_\{t\}is the multivariate input\. The LSTM unit uses forget \(ftf\_\{t\}\), input \(iti\_\{t\}\), and output \(oto\_\{t\}\) gates with candidate cell stateC~t\\tilde\{C\}\_\{t\}\. The GRU unit uses update \(ztz\_\{t\}\) and reset \(rtr\_\{t\}\) gates with candidate activationgtg\_\{t\}\. The LSTM carries forward long\- and short\-term memory separately as cell stateCCand hidden statehh, while GRU only forwards one combined state\. Blue and red frames indicate sigmoid \(σ\\sigma\) andtanh\\tanhactivations, respectively\. The recurrent weightsw,uw,uand biasesbbare learned\.The equations for the gated units and the forward pass of theLSTMandGRUare provided in Appendix[A\.1](https://arxiv.org/html/2610.19270#A1.SS1)\.

The general model architecture follows an RNN\-encoder design, specifically the encoder described by[Kirsch et al\. \[26\]](https://arxiv.org/html/2610.19270#bib.bib34)\. The encoder embeds the spectral\-temporal information into a higher\-dimensional hidden representation before progressively reducing it to retain the most relevant information\. We adapted this architecture for classification by replacing the decoder with a classification head, which maps the learned feature representation to the target classes\. In addition, we implemented an optionalMulti\-head Attention \(MHA\)\-based pooling mechanism, inspired by the Pooling by Multihead Attention \(PMA\) approach proposed by[Lee et al\. \[30\]](https://arxiv.org/html/2610.19270#bib.bib9), prior to the classification head\. Unlike classification approaches that rely solely on the final hidden state, the attention pooling mechanism incorporates information from all RNN time steps\. A learnable query attends to the temporal latent representations, assigning attention weights to the individual time steps and aggregating them into a single fixed\-dimensional representation that is subsequently passed to the classification head\. The architecture of an RNN model withMHApooling is illustrated in Figure[4](https://arxiv.org/html/2610.19270#S2.F4)\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/figs/GRU_MHA_Architekture_landsacpe_new3.png)Figure 4:GRU\-MHA Classifier Architecture\.The RNN output is projected into key and value representations, while a learnable global query provides the query representation\. After splitting the projected representations into multiple attention heads, scaled dot\-product attention is computed between the global query and the keys of all time steps\. The resulting attention scores are normalized using softmax and used to calculate a weighted sum of the value vectors, yielding one pooled representation per head\. These are connected and used as inputs for the classification head\.We incorporatedMHApooling to allow the model to weight the temporal representations produced by the recurrent encoder\. Because time steps contribute unequally to classification, theMHApooling enables the model to learn time steps containing particularly discriminative information\. By using multiple attention heads, different heads can learn complementary weighting patterns\[[31](https://arxiv.org/html/2610.19270#bib.bib23)\], potentially capturing information associated with distinct phenological periods\. We therefore expectedMHApooling to produce a more informative latent\-space representation without increasing the dimensionality of the latent representation\. Consequently, we expected it to improve classification, particularly for classes whose spectral differences are most pronounced during specific periods of the growing season\. This mechanism is schematically illustrated in Figure[5](https://arxiv.org/html/2610.19270#S2.F5)\. For simplicity, attention is shown for a single spectral band, although in the model it is applied to the full embedded representation of all bands at each time step\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/figs/attention_plot_Single_COL.png)Figure 5:Multi\-head attention weighting\.The top panel shows the temporal profile of spectral band 11 \(SWIR 1\), which is included for visualization purposes\. The panels below show the attention weights assigned to each time step by the individual attention heads\. The blue bars in the top panel represent the attention weights aggregated across all heads for each time step\. The attention mechanism operates on the latent representations rather than directly on the spectral band values\.Because the time series are incomplete and irregularly spaced, we included theday of year \(DOY\)as a scaled continuous input feature\. This provided the recurrent models with explicit calendar\-time information and relaxed the assumption that consecutive observations are equally spaced\. Because each sequence covered a single calendar year and did not span the December–January transition, cyclic encoding ofDOYwas not required\.

### 2\.3Feature Importance

Deep learning models are often considered black boxes because they provide limited insight into the factors driving their predictions\[[32](https://arxiv.org/html/2610.19270#bib.bib28),[33](https://arxiv.org/html/2610.19270#bib.bib29)\]\. There are several attribution methods that assign a scalar contribution value to each input feature, aiming to quantify how strongly each feature influenced the network’s prediction\[[34](https://arxiv.org/html/2610.19270#bib.bib24)\]\. These methods can make it difficult to determine whether observed inaccuracies originate from the model itself or from the attribution method\[[35](https://arxiv.org/html/2610.19270#bib.bib37)\]\. We selectedIntegrated Gradients \(IG\), introduced by[Sundararajan et al\. \[35\]](https://arxiv.org/html/2610.19270#bib.bib37), as the attribution method for this study\.IGis a widely used and theoretically grounded attribution method for explaining predictions of differentiable models\[[36](https://arxiv.org/html/2610.19270#bib.bib26),[37](https://arxiv.org/html/2610.19270#bib.bib25)\]\. It has also been shown to provide informative explanations for multivariate time\-series models\[[38](https://arxiv.org/html/2610.19270#bib.bib27)\]\.

Integrated Gradients is a gradient\-based attribution method that assigns an attribution value to each input feature, that is, to each spectral band at each time step\. For each input feature,IGintegrates the gradients of the model output along a linear interpolation path between a predefined baseline and the observed multivariate time series\[[35](https://arxiv.org/html/2610.19270#bib.bib37)\]\. The resulting attribution values quantify how each feature contributes to the difference between the model output for the observed input and that for the baseline\. Consequently, the resulting attribution matrix represents the contribution of each spectral\-band\-time\-step combination to the prediction of the respective target class relative to the baseline input\.

IGi\(x\)::=\(xi−xi′\)×∫α=01=∂F⁡\(x′\+α×\(x−x′\)\)∂xidα\\operatorname\{IG\}\_\{i\}\(x\)::=\(x\_\{i\}\-x^\{\\prime\}\_\{i\}\)\\times\\int^\{1\}\_\{\\alpha=0\}=\\frac\{\\partial F\(x^\{\\prime\}\+\\alpha\\times\(x\-x^\{\\prime\}\)\)\}\{\\partial x\_\{i\}\}d\\alpha\(1\)
Herexxdenotes the input sample andx′x^\{\\prime\}denotes the corresponding baseline input\. The term∂F⁡\(x\)∂xi\\frac\{\\partial F\(x\)\}\{\\partial x\_\{i\}\}is the gradient of the model outputF⁡\(x\)F\(x\)along theit​hi^\{th\}dimension\. In practice, the integral is approximated using a Riemann sum withmmsteps\.

IGiapprox⁡\(x\)::=\(xi−xi′\)×∑k=1m∂F⁡\(x′\+km×\(x−x′\)\)∂xi×1m\\operatorname\{IG\}^\{\\mathrm\{approx\}\}\_\{i\}\(x\)::=\(x\_\{i\}\-x^\{\\prime\}\_\{i\}\)\\times\\sum\_\{k=1\}^\{m\}\\frac\{\\partial F\\left\(x^\{\\prime\}\+\\frac\{k\}\{m\}\\times\(x\-x^\{\\prime\}\)\\right\)\}\{\\partial x\_\{i\}\}\\times\\frac\{1\}\{m\}\(2\)
The authors suggest comparing observations with a zero baseline to identify the input features most relevant to a model’s classification decision\[[35](https://arxiv.org/html/2610.19270#bib.bib37)\]\. By applying IG to multiple observations and averaging the results, a more robust estimate of feature importance can be obtained\.

Because we focused on classes that were difficult to distinguish, we also used representative baselines derived from the respective comparison classes\. These class\-specific baselines allowed us to identify the spectral bands and time steps that contributed most strongly to distinguishing the observed class from the class used as the baseline\.

## 3Methodology

In this section, we describe the methodological framework used for forest classification and model interpretation\. We first describe the data preprocessing steps, then the training and inference procedures for the recurrent neural network and ensemble models\. We then introduce the evaluation procedure and the Integrated Gradients approach used to analyze the spectral and temporal relevance of the model inputs\. To support reproducibility, the code used in this study for data preprocessing, model development, evaluation, and Integrated Gradients analysis is publicly available on GitHub111GitHub repository:[https://github\.com/Max25832/TreeClassXAI](https://github.com/Max25832/TreeClassXAI)\.

### 3\.1Data Preprocessing

Many recent studies use multi\-year satellite image time series for forest classification\[[14](https://arxiv.org/html/2610.19270#bib.bib33),[10](https://arxiv.org/html/2610.19270#bib.bib19),[18](https://arxiv.org/html/2610.19270#bib.bib14)\]\. In contrast, we focused on single\-season sequences because the objective of this study is to identify phenologically relevant time periods and spectral regions for separating similar forest structures, rather than to optimize an operational tree species mapping model\.

The time series were therefore divided into yearly sequences that can be generalized across annual phenological developments\. This reduced the risk of the models learning patterns specific to individual extreme events, such as the exceptionally hot and dry summer of 2018\. Furthermore, this approach also increased the number of available training examples\. Additionally, we excluded observations before 2017, because Sentinel\-2B started in June 2017, leading to more observations per year\[[20](https://arxiv.org/html/2610.19270#bib.bib36)\]\. The train\-test split was performed prior to temporal segmentation to avoid data leakage\. This ensured that time series segments of the same pixel from different years were strictly assigned to either the training or the test dataset and did not appear in both datasets simultaneously\. We used a standard 80:20 train\-test split\. To reduce class imbalance in the training data, the number of samples per class was limited to at most twice the sample size of the smallest class:

nc≤2⋅mink∈𝒞⁡nk∀c∈𝒞n\_\{c\}\\leq 2\\cdot\\min\_\{k\\in\\mathcal\{C\}\}n\_\{k\}\\quad\\forall c\\in\\mathcal\{C\}\(3\)wherencn\_\{c\}denotes the number of samples in classcc,𝒞\\mathcal\{C\}denotes the set of all classes,kkis an index over the classes in𝒞\\mathcal\{C\}, andmink∈𝒞⁡nk\\min\_\{k\\in\\mathcal\{C\}\}n\_\{k\}denotes the sample size of the smallest class\. Downsampling was performed after the data were divided into yearly subsequences\. This allowed us to retain a spatially and temporally diverse subset of the highly represented classes, as samples were randomly selected from different pixels and years, thereby increasing both spatial and temporal diversity\.

Table 1:Final Sequences Sample SizeBecause the deep learning models used sequential inputs, the original time\-series structure could be retained\. In contrast, the ensemble models required tabular input data\. For these models, the time series were transformed into a wide\-format representation, where each spectral band at each time step was represented by a separate predictor variable\. Using daily time steps would have resulted in366​days×10​spectral bands=3,660366\\text\{ days\}\\times 10\\text\{ spectral bands\}=3\{,\}660predictors\. This would have produced a very high\-dimensional dataset, increasing the risk of overfitting and leading to a large number of missing values, as the time series contained, on average, only approximately 24 observations per year\. To reduce dimensionality and obtain a more compact representation, we aggregated the data to monthly median values\. Remaining gaps were imputed by forward filling\. As a result, the tabular dataset used for the ensemble models contained 120 predictors\.

### 3\.2Model Training and Inference

Hyperparameter tuning for both ensemble models was performed on the training set using five\-fold cross\-validation and predefined tuning grids\. The evaluated hyperparameters are provided in Appendix Table[6](https://arxiv.org/html/2610.19270#A1.T6)\. For the deep learning models, the training set was further split into training and internal validation subsets using an 80:20 split\. The internal validation subset was used for checkpoint selection during training, while the hold\-out test set was excluded from training, validation, and hyperparameter tuning\. After evaluation, test\-set Macro\-F1 was used to select the model configurations for the feature importance analysis described below\. All RNN\-based models were trained for 80 epochs using the Adam optimizer with a learning rate of1×10−31\\times 10^\{\-3\}\. For the loss function, we used a weighted combination of cross\-entropy loss and a Macro\-F1\-based loss, with weights of 0\.7 for cross\-entropy loss and 0\.3 for the Macro\-F1\-based loss\. We additionally applied early stopping during training, using Macro\-F1 on the validation set as the monitoring metric\. Both recurrent architectures, GRU and LSTM, were trained with and without multi\-head attention pooling\. All RNN\-based models used a bidirectional variant, in which the input sequence was processed in both forward and backward temporal directions\. We tested bottleneck dimensions of 32, 64, and 128\. The initial embedding dimension was set to four times the corresponding bottleneck dimension \(see Appendix Table[6](https://arxiv.org/html/2610.19270#A1.T6)\)\.

After training, all models were evaluated on the held\-out test set using accuracy, balanced accuracy, and Macro\-F1\. Balanced accuracy and Macro\-F1 were included to account for class imbalance in the dataset\.

### 3\.3Application to a Mixed\-Forest Area

The dataset authors note that it is intended primarily for model training and is less suitable for independent validation\[[20](https://arxiv.org/html/2610.19270#bib.bib36)\]\. To account for this and to illustrate the practical application of the trained model, we visually selected an additional mixed\-forest area in Thuringia from orthorectified aerial images222Orthophotos:[https://geoportal\.thueringen\.de/gdi\-th/download\-offene\-geodaten/download\-luftbilder\-und\-orthophotos](https://geoportal.thueringen.de/gdi-th/download-offene-geodaten/download-luftbilder-und-orthophotos)\.\. We extracted Sentinel\-2 time series for each pixel within the area using Google Earth Engine and classified them with the best\-performing trained model\. The resulting pixel\-wise predictions were used to generate a species\-composition map and to estimate the area assigned to each forest class\. We used orthophotos from 2023 for visual confirmation and classified the years 2022 and 2023\. For each pixel, we selected the year with the higher classification confidence to generate the final map\.

### 3\.4Feature Importance Analysis

To investigate model explainability, we appliedIGto quantify the contribution of individual spectral bands and time steps to each classification\. For each of the main tree species, we initially computed attributions relative to a zero baseline\. We randomly selected 2,000 multivariate time series per class from the test set and restricted the analysis to samples that were correctly classified by the respective model\. This was done to ensure that the resulting attribution patterns reflected successful model decisions\. The path integral inIGwas approximated using a Riemann sum withm=30m=30interpolation steps between the baseline and the observed input\. Furthermore, we compared observed time series against class\-specific baselines derived from other classes to identify the spectral and temporal features that contributed most strongly to distinguishing the observed class from the respective baseline class\. In this context, an individual input feature corresponds to the value of a single spectral band at a single time step\. BecauseIGrequires the baseline and input to have the same shape, we had to account for the varying sequence lengths\. We therefore constructed a median baseline time series for each selected baseline species by calculating the median separately for each spectral band and day of year across all available observations\. For each observed sequence, we then selected only those values from the median baseline that corresponded to the day of year present in that sequence, thereby ensuring temporal alignment between the input and its baseline\. This process is illustrated for a single spectral band of one observation in Figure[6](https://arxiv.org/html/2610.19270#S3.F6)\.

1000200030004000B8Baseline for the Integrated Gradeints Class ComparisonMedian Baseline Species020004000B8Observation Timeseries050100150200250300350Day of Year01000200030004000B8New Baseline

Figure 6:Class\-specific baseline derivation forIGattribution\.For each observation, a baseline with the same temporal structure was derived from the median time series of the comparison class\. The example shows an oak observation \(blue\) and the corresponding beech baseline \(green\) for band 11, derived from the median of all beech observations \(red\) and matched to the observation time steps\.For the feature importance analysis, we selected the model configurations based on the highest test\-set Macro\-F1 scores\. We did not compare models that differed only in embedding dimensions\. In contrast, we separately analyzed models with and without multi\-head attention pooling because these architectures aggregate the temporal representations differently and may therefore differ in attribution patterns\. We compared all classes against the zero baseline\. In addition, pairwise class\-specific baseline comparisons were performed betweenpure Norway spruceandmixed Norway spruceand betweenoakandbeech\. Each pairwise comparison was conducted in both directions, with each class serving once as the observed class and once as the baseline class\.

## 4Results

To evaluate model performance, we comparedRNNconfigurations with different recurrent architectures, input and bottleneck dimensions, and eitherMHApooling or no pooling against the ensemble models\. Performance was assessed primarily using the Macro\-F1 score, supported by additional classification metrics\. We further appliedIGto analyze the spectral and temporal patterns contributing to correct class predictions and conducted detailed class\-comparison analyses for Norway spruce versus mixed Norway spruce and oak versus beech\.

### 4\.1Classification Performance

Based on the Macro\-F1 score, the GRU model with multi\-head attention and a bottleneck dimension of 128 achieved the best overall performance, with a Macro\-F1 score of 0\.851\. It additionally reached an accuracy of 0\.854 and a balanced accuracy of 0\.933\. The GRU\-128 model without attention performed similarly but achieved a slightly lower Macro\-F1 score of 0\.847, together with an accuracy of 0\.853 and the same balanced accuracy of 0\.933\. Among the baseline models, XGBoost achieved the highest Macro\-F1 score of 0\.839 and an accuracy of 0\.851, whereas the Random Forest showed the lowest overall performance\. For both recurrent architectures, larger bottleneck dimensions generally corresponded to higher Macro\-F1 scores, while the effect of multi\-head attention was comparatively small\.MHAdid not improve the Macro\-F1 performance of theLSTMmodels, but increased Macro\-F1 by 0\.004 for theGRUmodels\. These results are shown in Table[2](https://arxiv.org/html/2610.19270#S4.T2)\.

Table 2:Comparison of Model Configurations\* BN Dim\. = Model Bottleneck Dimension; Acc = accuracy; BA = balanced accuracy\.

The confusion matrices in Figure[7](https://arxiv.org/html/2610.19270#S4.F7)revealed that the main classification challenge was distinguishing Norway spruce from mixed Norway spruce\. Across all models, approximately one\-fifth of the samples in each class were confused with the other class\. The green confusion matrices are normalized by true class and report the proportion of samples in percent; therefore, the diagonal entries correspond to class\-specific sensitivity\. For the GRU\-MHA model, 21\.6% of Norway spruce samples were classified as mixed Norway spruce, while 21\.2% of mixed Norway spruce samples were assigned to Norway spruce\. Accordingly, class\-specific sensitivity was comparatively low for Norway spruce \(76\.3%\) and mixed Norway spruce \(72\.3%\)\. The remaining classes achieved substantially higher sensitivities\. The GRU\-MHA model reached 97\.2% for Scots pine, 96\.3% for beech, 90\.5% for disturbed Norway spruce, 93\.7% for oak, and 97\.8% for soil\. Compared with the corresponding GRU model without MHA, the addition of MHA reduced the proportion of oak samples misclassified as beech from 4\.5% to 2\.6%, while the reverse confusion increased slightly from 1\.3% to 1\.8%\. This indicates that MHA produced only small, class\-specific changes in the confusion patterns\. The mutual confusion between pure and mixed Norway spruce remained the dominant source of classification error\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/confusionmatrices_sc_percentage-img0.png)NWSNWS MixedScots PineBeechNWS DisturbedOakSoil14707455511718201147681392949478611913594131299552731417936642284674664112919243970920811623413101720582941364153899617636GRU![Refer to caption](https://arxiv.org/html/2610.19270v1/confusionmatrices_sc_percentage-img1.png)1502542501532017217518581292819449212714011129738531584630472364107365761012083645109102018394397593210329754542605910817603GRU\-MHA![Refer to caption](https://arxiv.org/html/2610.19270v1/confusionmatrices_sc_percentage-img2.png)118713296143142609577443238087909323203089457476283033192262509845648158685065179131748818029101055177831162915491008115019XGBoost![Refer to caption](https://arxiv.org/html/2610.19270v1/confusionmatrices_sc_percentage-img3.png)NWSNWS MixedScots PineBeechNWS DisturbedOakSoilNWSNWS MixedScots PineBeechNWS DisturbedOakSoil74\.7%23\.1%0\.6%0\.1%1\.0%0\.1%0\.4%20\.1%72\.8%1\.9%2\.9%0\.9%1\.0%0\.3%0\.3%1\.6%96\.7%0\.0%1\.1%0\.1%0\.1%0\.1%1\.0%0\.1%97\.3%0\.2%1\.3%0\.1%1\.1%1\.7%3\.1%0\.4%93\.3%0\.0%0\.3%0\.1%1\.5%0\.6%4\.5%0\.3%91\.7%1\.3%0\.2%0\.2%0\.2%0\.3%0\.5%0\.5%98\.0%![Refer to caption](https://arxiv.org/html/2610.19270v1/confusionmatrices_sc_percentage-img4.png)NWSNWS MixedScots PineBeechNWS DisturbedOakSoil76\.3%21\.6%0\.8%0\.1%0\.9%0\.1%0\.3%21\.2%72\.3%2\.3%2\.3%0\.7%1\.0%0\.3%0\.3%1\.2%97\.2%0\.0%0\.9%0\.2%0\.1%0\.1%1\.6%0\.0%96\.3%0\.1%1\.8%0\.1%1\.6%2\.0%4\.9%0\.4%90\.5%0\.1%0\.4%0\.2%1\.7%0\.3%2\.6%0\.1%93\.7%1\.3%0\.4%0\.3%0\.2%0\.3%0\.3%0\.6%97\.8%![Refer to caption](https://arxiv.org/html/2610.19270v1/confusionmatrices_sc_percentage-img5.png)NWSNWS MixedScots PineBeechNWS DisturbedOakSoil75\.9%21\.1%0\.9%0\.1%1\.7%0\.1%0\.4%22\.1%70\.7%2\.3%2\.8%0\.9%0\.9%0\.3%0\.2%1\.6%97\.3%0\.0%0\.7%0\.1%0\.1%0\.0%1\.7%0\.1%96\.4%0\.3%1\.5%0\.1%2\.4%3\.1%8\.6%0\.6%84\.0%0\.4%0\.9%0\.0%1\.5%0\.5%5\.4%0\.3%90\.8%1\.6%0\.1%0\.2%0\.1%0\.3%0\.7%0\.5%98\.1%Predicted classTrue class

Figure 7:Confusion matrices for the two best\-performing RNN models, GRU with and without attention, and the XGBoost modelThe class\-specific evaluation metrics support the patterns observed in the confusion matrices\. Norway spruce achieved the lowest F1 scores, ranging from 0\.614 to 0\.691\. The GRU\-128 model with multi\-head attention achieved the highest F1 score \(0\.691\) and sensitivity \(0\.763\), while the highest specificity was 0\.919\. For mixed Norway spruce, F1 scores ranged from 0\.697 to 0\.784, with two recurrent models tying for the highest score\. The maximum sensitivity and specificity were 0\.729 and 0\.943, respectively\. For both classes, specificity exceeded sensitivity, indicating that false\-negative rates were higher than false\-positive rates\. The remaining classes achieved higher F1 scores\. The best scores were obtained by the GRU\-128 model with multi\-head attention for disturbed areas \(0\.797\), the GRU\-128 model without attention for Scots pine \(0\.968\), and the LSTM\-128 model without attention for beech \(0\.909\)\. XGBoost achieved the highest scores for oak \(0\.836\) and soil \(0\.983\), although the recurrent models performed similarly\. Overall, soil and Scots pine had the highest class\-specific F1 scores, whereas Norway spruce had the lowest\. Random Forest generally underperformed compared with XGBoost and the recurrent models\. These results are presented in detail in Table[3](https://arxiv.org/html/2610.19270#S4.T3)\.

Table 3:Per\-class comparison of model configurations
### 4\.2Performance in a Mixed\-Forest Area

To additionally assess model performance in a mixed\-forest area, we applied the GRU\-MHA\-128 model to single\-year time series from a randomly selected site in Thuringia\. The resulting classification map is shown in Fig\.[8](https://arxiv.org/html/2610.19270#S4.F8)\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/figs/ClassificationMap.png)\(a\)Predicted forest classes
11\.48 ha0\.05 ha2\.32 ha2\.395 ha0\.01 ha03006009001200ClassPixel CountClassBeechNorway Spruce: DisturbedNorway Spruce: HealthyNorway Spruce: Mixedsoil

\(b\)Classification statistics

Figure 8:Map with distribution statistics of predicted forest classes\.The model classified 11\.48 ha as beech, 2\.32 ha as pure Norway spruce, 2\.40 ha as mixed Norway spruce, 0\.05 ha as disturbed forest, and 0\.01 ha as open soil\. Because pixel\-level labels were not available, we used an orthorectified aerial image as a qualitative visual reference\. The visual comparison indicated reasonable agreement between the predicted classes and the aerial image, providing an additional qualitative assessment of model performance outside the original dataset\.

### 4\.3Model Explainability

To improve the explainability of the deep learning model, we appliedIGin two ways\. First, individual class attributions were calculated using a zero baseline\. Second, we directly compared Norway spruce with mixed Norway spruce and oak with beech by using the respective opposing class as the baseline\. The results of these analyses are presented in the following sections\.

#### 4\.3\.1Class\-Specific Integrated Gradients

Hereafter, spectral bands are referred to by their band identifiers, as defined in Table[4](https://arxiv.org/html/2610.19270#S4.T4)\.

Table 4:Sentinel\-2 Level\-2A bands used in this study\.\* NIR = near infrared; SWIR = short\-wave infrared\.

The feature attribution analysis revealed distinct temporal and spectral patterns among the classes\. For beech, attribution values peaked aroundDOY120 in bands B3, B5, B6, and B11\. Oak showed a similar temporal pattern, but with lower attribution values in bands B3–B6 and an earlier peak in B12, approximately betweenDOY110 and 125\. For Scots pine, theIGanalysis showed high B11 attribution values between approximatelyDOY130 and 190, whereas B6 exhibited moderate attribution values throughout the year\. Pure Norway spruce and mixed Norway spruce showed similar attribution patterns, with high values in B6 aroundDOY150\. For mixed Norway spruce, the B6 attributions were concentrated between approximatelyDOY140 and 175\. In contrast, for pure Norway spruce, they extended from approximatelyDOY120 to 220\. Both Norway spruce and Scots pine showed elevated attribution values between approximatelyDOY100 and 180\. However, Scots pine also exhibited elevated attributions toward the end of the year\. The disturbed Norway spruce class showed no pronounced temporal attribution pattern\. Attribution values were broadly similar across bands, except for B2, B7, and B12\.

Overall, theIGanalysis revealed consistent differences in attribution magnitude among the spectral bands\. Across classes, B2 and B12 generally showed the lowest attribution values, whereas B6 and B11 showed the highest attributions for correct species classifications\. The results are visualized in Figure[9](https://arxiv.org/html/2610.19270#S4.F9)\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_0BL-img0.png)50100150200250300350DOY02Mean \|IG\|Scots Pine0\.1001\.1812\.262Mean \|IG\| / BandB2B3B4B5B6B7B8B8AB11B12![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_0BL-img1.png)50100150200250300350DOY02Mean \|IG\|Beech0\.1001\.1812\.262Mean \|IG\| / BandB2B3B4B5B6B7B8B8AB11B12![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_0BL-img2.png)50100150200250300350DOY02Mean \|IG\|Oak0\.1001\.1812\.262Mean \|IG\| / BandB2B3B4B5B6B7B8B8AB11B12![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_0BL-img3.png)50100150200250300350DOY02Mean \|IG\|Norway Spruce0\.1001\.1812\.262Mean \|IG\| / BandB2B3B4B5B6B7B8B8AB11B12![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_0BL-img4.png)50100150200250300350DOY02Mean \|IG\|Norway Spruce mixed0\.1001\.1812\.262Mean \|IG\| / BandB2B3B4B5B6B7B8B8AB11B12![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_0BL-img5.png)50100150200250300350DOY02Mean \|IG\|Norway Spruce Disturbed0\.1001\.1812\.262Mean \|IG\| / BandB2B3B4B5B6B7B8B8AB11B12![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_0BL-img6.png)01234567Mean absolute attribution

Figure 9:Integrated Gradients attribution results for the six investigated classes against a zero\-baseline\.
#### 4\.3\.2Integrated Gradients for Class Comparison

TheIGanalysis focused on the distinction between pure and mixed Norway spruce stands and between beech and oak\. Because beech and oak were classified with comparatively high performance, the following analysis focuses on the attribution patterns for pure and mixed Norway spruce, while also examining differences between beech and oak\.

In the comparison between pure and mixed Norway spruce, B6 and B11 showed the highest overall attribution values\. Most other bands exhibited comparable cumulative attributions across the year, whereas B2 and B12 consistently showed lower values\. B11 exhibited elevated attributions throughout much of the year, with a pronounced peak between approximatelyDOY140 and 180\. Attributions for B6 increased after approximatelyDOY125 and remained elevated until aroundDOY290\. The visible bands B2, B3, and B4 and the red\-edge band B5 showed comparatively low but persistent attributions throughout the year\. Bands B7, B8, and B8A followed a temporal pattern similar to that of B6, although their attribution magnitudes were generally lower\. Overall, the highest attribution values for distinguishing pure from mixed Norway spruce occurred between approximatelyDOY125 and 225\. This comparison is visualized in Figure[10](https://arxiv.org/html/2610.19270#S4.F10)\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_NWS_mixedBL_1Col-img0.png)50100150200250300350DOY0\.00\.5Mean \|IG\|Model Attributions for Pure Norway Spruce using Mixed as Baseline0\.00\.51\.0Mean \|IG\| / BandB2B3B4B5B6B7B8B8AB11B12![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_NWS_mixedBL_1Col-img1.png)0\.000\.250\.500\.751\.001\.251\.501\.752\.00Mean absolute attribution

Figure 10:Attribution map for the separation of pure and mixed Norway spruce\.For beech versus oak, B3, B5, and B6 showed attribution patterns similar to those observed for beech relative to the zero baseline\. However, B6 exhibited elevated attributions over a broader period, approximately fromDOY110 to 140\. The highest attribution values were observed for B6 and B11, particularly between approximatelyDOY110 and 150\. This comparison is visualized in Figure[11](https://arxiv.org/html/2610.19270#S4.F11)\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_Beech_oakBL_1Col-img0.png)50100150200250300350DOY0\.00\.5Mean \|IG\|Model Attributions for Beech using Oak as Baseline0\.00\.250\.5Mean \|IG\| / BandB2B3B4B5B6B7B8B8AB11B12![Refer to caption](https://arxiv.org/html/2610.19270v1/attributions_Beech_oakBL_1Col-img1.png)0\.00\.20\.40\.60\.81\.0Mean absolute attribution

Figure 11:Attribution Map for the separation of beech and oak

## 5Discussion

This study demonstrates the value of recurrent neural networks and attribution\-based analysis for classifying and interpreting annual multispectral time series\. The RNN models outperformed XGBoost, highlighting their ability to capture temporal dependencies in the data\. Combined with Integrated Gradients, they also enabled a more fine\-grained interpretation of relevant spectral bands and time periods than feature\-importance measures from tree\-based models, which are typically aggregated over broader temporal intervals\. In addition, the class\-comparison attribution analysis provided insights beyond the conventional zero\-baseline approach by revealing shifts in the relative importance of individual spectral bands\. These findings show that directly comparing classes can improve the interpretation of spectral and temporal differences between closely related forest classes\.

### 5\.1Comparative Model Performance

The class\-specific results indicate that classification performance was primarily determined by the separability of the spectral\-temporal classes rather than by the architecture of the model\. The main source of error across all configurations was the mutual confusion between pure and mixed Norway spruce, with approximately one\-fifth of the observations misclassified between these two classes\. This persistent confusion is likely due to the strong overlap in multispectral and temporal characteristics between the two classes\. Pure stands tend to occupy more distinct regions of spectral\-temporal feature space, whereas mixed pixels are strongly influenced by their species composition and form a gradual spectral continuum\. The difficulty of separating these classes therefore appears to be constrained not only by model capacity but also by the information content of the Sentinel\-2 time series and the definition of the classes themselves\. In contrast, Scots pine, beech, and soil were consistently classified with high accuracy across model architectures, further supporting the interpretation that class separability was the dominant factor determining classification success\. Disturbed Norway spruce was also separated more clearly from healthy spruce than pure and mixed healthy spruce were from each other\. Disturbance alters canopy condition and therefore produces stronger deviations in spectral\-temporal trajectories, which may produce a more distinctive spectral\-temporal signal than differences in species composition within otherwise healthy spruce stands\.

Despite these class\-dependent limitations, the recurrent models achieved the highest overall classification performance, with the GRU using multi\-head attention and a bottleneck dimension of 128 reaching the highest Macro\-F1 score of 0\.851, while Random Forest consistently performed worst\. However, the improvement over XGBoost was relatively small\. This suggests that, although explicit sequence modeling provides an advantage for annual multispectral time\-series classification, much of the relevant temporal information can also be captured through engineered or temporally aggregated features\. The benefit of recurrent modeling was particularly evident for disturbed Norway spruce, where the recurrent models achieved substantially higher F1 scores and sensitivity than XGBoost, indicating that explicit modeling of temporal dependencies may be especially useful for classes characterized by changes in canopy condition over time\. In contrast, improvements for most healthy tree\-species classes were small, and XGBoost remained competitive or even performed best for some well\-separated classes\.

Larger bottleneck dimensions generally improved the performance of both recurrent architectures, suggesting that sufficient latent capacity may be more important than the choice between LSTM and GRU or the inclusion of attention\. Although larger bottleneck dimensions may improve classification by preserving more information, they also increase the risk of overfitting\. This trade\-off should be considered when selecting the latent representation size\. The funnel\-shaped encoder progressively compresses the hidden representation and may therefore encourage the model to retain the most relevant information\. However, the present experiments cannot determine whether the funnel architecture itself is beneficial or whether the observed gains primarily result from the increased latent dimensionality and overall model capacity\. A direct comparison with a non\-funneling architecture of similar size should be considered in future work\.

MHAhad only a limited effect on predictive performance\. It produced a marginal improvement for the GRU\-128 model and did not improve the Macro\-F1 score of the LSTM\-128 model, suggesting that the recurrent hidden states already encoded most of the temporal information required for classification\. Nevertheless,MHApooling improved the Macro\-F1 score across all GRU configurations, with somewhat larger gains for models with smaller bottleneck dimensions, and produced minor class\-specific improvements\. Attention may therefore compensate for limited latent capacity by suppressing less informative temporal representations\. However, its main value in this study lies in improving the temporal interpretation of the model by emphasizing relevant time steps and making their contribution more apparent in the attribution analysis\.

Overall, the performance gains over XGBoost were class\-dependent and moderate\. Latent\-space size had a clearer influence on performance thanMHA, while the remaining errors were driven mainly by limited spectral\-temporal separability rather than by the selected model architecture\. The main advantage of the RNNs, however, lies in their temporally structured representation, which supports a more detailed interpretation of phenologically relevant periods\.

### 5\.2Explainability

The class\-comparisonIGanalysis provided insights beyond the zero\-baseline analysis\. Relative to a zero baseline, B6 showed the highest attributions for both pure and mixed Norway spruce\. However, when each class was evaluated against a baseline derived from the other class, B11 became the most strongly attributed band\. Thus, the features supporting individual class predictions differed from those used to distinguish the two closely related classes\. The direct comparison also revealed a seasonal pattern that was not apparent in the zero\-baseline results\. Attributions for B6, B7, B8 and B8A increased after approximatelyDOY125, broadly coinciding with leaf development in Central European broadleaf species\[[39](https://arxiv.org/html/2610.19270#bib.bib7)\]\. The differences become clearer when the median reflectance time series of pure and mixed Norway spruce are compared \(see Figure[12](https://arxiv.org/html/2610.19270#S5.F12)\)\. The median B6 reflectance trajectories of the two classes begin to diverge aroundDOY100 and remain separated for most of the growing season\. However, the highest attribution values during this period occur in B11 rather than B6\. In particular, B11 attributions increased when the difference between the B6 trajectories was most pronounced\. B11 reflectance also remained consistently higher for the mixed class throughout the year, with a peak in late spring\. These results suggest that the model did not rely on individual bands independently but rather on their combined spectral\-temporal patterns when distinguishing pure from mixed Norway spruce stands\.

![Refer to caption](https://arxiv.org/html/2610.19270v1/figs/nws_mixed_b6_b11_MHA.png)Figure 12:Separation between pure and mixed Norway spruce\.The top row shows the median test\-set time series with confidence bands\. The bottom row displays Integrated Gradients attributions along the Norway spruce median time series, highlighting time steps relevant for distinguishing the two species\.This result is particularly relevant because pure and mixed Norway spruce showed the highest degree of mutual misclassification across the evaluated models\. Although the generally high relevance of red\-edge and short\-wave infrared bands is consistent with previous remote\-sensing research\[[14](https://arxiv.org/html/2610.19270#bib.bib33),[22](https://arxiv.org/html/2610.19270#bib.bib18)\], the class\-comparison temporal attribution patterns provide more detailed insight into how the models distinguished these closely related forest classes\.

Although the phenological patterns of beech and oak appeared more similar than those of pure and mixed Norway spruce, the models distinguished the two broadleaf species with high accuracy\. However, the attribution patterns did not correspond directly to the periods with the largest visible differences between the median time series\. Instead, the highest attributions occurred mainly betweenDOY100 and 150, coinciding with the leaf\-unfolding period\. Although the median time series showed little separation between beech and oak during this period, the model appears to have used subtle spectral–temporal differences to distinguish the two species \(Figure[13](https://arxiv.org/html/2610.19270#S5.F13)\)\.

020040060080010001200B31000200030004000B65001000150020002500B11100200300200400600800100012001002003001000200030001002003005001000150020002500![Refer to caption](https://arxiv.org/html/2610.19270v1/oak_beech_b3_b6_b11_MHA-img0.png)0\.250\.500\.751\.00Attribution strengthDOYReflectanceOak MedianOak 2nd–98th percentileBeech MedianBeech 2nd–98th percentile

Figure 13:Seperation between Beech and Oak\. The top row shows the median test\-set time series with confidence bands\. The bottom row displays Integrated Gradients attributions along the beech median time series, highlighting time steps relevant for distinguishing the two species\.
### 5\.3Potential Applications in Forestry

We expected the spectral\-temporal patterns of species to vary severely between years, potentially reducing the effectiveness of single\-year time series\. In particular, extreme summers such as 2018 and 2022 in Germany were expected to alter spectral patterns and increase misclassification\. However, the proportion of false predictions did not increase noticeably in these years \(Fig\.[14](https://arxiv.org/html/2610.19270#S5.F14)\)\.

BeechOakPineNWSNWS mixedNWS dist\.Soil2017201820192020202120222017201820192020202120222017201820192020202120222017201820192020202120222017201820192020202120222017201820192020202120222017201820192020202120220%25%50%75%100%YearPercentagePredictionCorrectFalseError bars: ± 1 SD across models

Figure 14:False Predictions over years:Here the mean percentage±\\pmSD of false predictions for each class is visualized for every observed year across all RNN models\. Norway spruce is abbreviated toNWS\. The disturbed Norway spruce class has no observations for the years 2021 and 2022, as there was no data where the disturbance year was after 2020\.Across all RNN models, fewer annual observations generally reduced classification performance\. The impact of observation count on correct classification percentage for pure and mixed Norway spruce over all RNN models is visualized in Figure[15](https://arxiv.org/html/2610.19270#S5.F15)\. This offers potential for future studues, as the launch of Sentinel\-2C increased yearly observations in europe by about 33% and reduced the average datagap lenght\[[40](https://arxiv.org/html/2610.19270#bib.bib1)\], which suggests better interpretability with our approach\.

False Predictions010203040500%25%50%75%100%ObservationsPercentagePredictionCorrectFalse

Figure 15:Mean percentage of false predictions \(±\\pmSD\) for pure and mixed Norway spruce combined, grouped by the number of time steps per time series\.A direct comparison with studies using multi\-year time series for similar forest classes in Germany is difficult, as these studies often included a larger number of classes\[[14](https://arxiv.org/html/2610.19270#bib.bib33),[13](https://arxiv.org/html/2610.19270#bib.bib17),[18](https://arxiv.org/html/2610.19270#bib.bib14)\]\. This increases the potential for misclassification, particularly among broadleaf species such as beech and oak\. At the same time, other studies did not distinguish between pure and mixed Norway spruce stands, which likely contributed to higher Macro\-F1 scores for Norway spruce\. The classification tasks therefore differ substantially, and a robust comparison would require further targeted studies\. Nevertheless, the strong classification performance achieved here indicates that the appropriate input length is task\-specific and that single\-year time series should not be excluded by default\. For detailed forest mapping, longer multi\-year time series may provide higher classification performance\[[14](https://arxiv.org/html/2610.19270#bib.bib33),[13](https://arxiv.org/html/2610.19270#bib.bib17),[18](https://arxiv.org/html/2610.19270#bib.bib14)\]\. In contrast, shorter time series may be preferable when lightweight operational models or energy\-efficient training and inference are priorities\. Previous studies have shown that reducing the number of input observations can improve computational and energy efficiency\[[41](https://arxiv.org/html/2610.19270#bib.bib6)\]\.In such a setting, reducing the number of input features also becomes important, as it further contributes to better energy efficiency\[[41](https://arxiv.org/html/2610.19270#bib.bib6)\]\. TheIGanalysis indicated that bands B2 and B12 contributed little to the model predictions and could therefore be tested for removal to improve efficiency\. Reducing the input dimensionality may also allow similar performance with smaller hidden representations and fewer trainable parameters, thereby lowering computational demand and associated emissions\. However, this would require further evaluation\. Identifying the temporal and spectral observations that contribute most strongly to distinguishing pure and mixed stands provides useful information for forest composition mapping\. The insights of the explainability analysis could guide the development of future models that better separate pure and mixed stands by focusing on the most informative spectral bands and seasonal periods\. This may also enable feature and time\-step reduction, supporting the development of more lightweight and computationally efficient models\.

### 5\.4Limitations

The training dataset included observations from different locations and years, allowing the models to capture spatial and interannual variability\. However, this variability was limited to the period from 2017 to 2022\. In addition, the dataset providers noted that the dataset is well suited for model training but less suitable for independent validation\[[20](https://arxiv.org/html/2610.19270#bib.bib36)\]\. The reported classification performance should therefore be interpreted with caution\. Furhtermore, Spatial and temporal autocorrelation may have influenced the performance estimates despite measures taken to prevent data leakage\. To avoid any data leakage a leave\-one\-region\-out cross\-validation approach would be most appropriate\. In this setup the data is spatially divided into regions and the model is trained on all but one region\. The left out region is used for testing\. This is repeated iteratively until each region has served once as the test set\. We partly addressed these limitations by applying the model to an independent mixed\-forest region\. However, this provided only a limited visual assessment of model performance and cannot replace evaluation on a larger and more diverse independent test dataset\. Another limitation is that the study region is limited to Germany and the results may not transfer directly to other forest types, climatic conditions or management systems\. We selected species classes based on the four predominant tree species in Germany, with different health stages and structures for Norway spruce, which also sets us apart from other german studies in terms of our class definitions and study design\. This makes it difficult to draw direct comparisons with previous studies regarding model performance, particularly with regard to the separability of pure and mixed spruce stands, as other studies have not focused on this approach\. Finally, the potential efficiency gains from shorter time series, feature reduction, and smaller model architectures were not measured directly\. Their effects on computational demand, energy consumption, and classification performance therefore require further evaluation\.

## 6Conclusion

This study investigated the classification and interpretation of annual Sentinel\-2 time series for forest species mapping, with particular emphasis on the discrimination of spectrally and phenologically similar classes\. Recurrent deep learning models were compared with established tree\-based ensemble methods, while their spectral\-temporal decision processes were examined using Integrated Gradients\. The best\-performing model, GRU\-MHA\-128, achieved a Macro\-F1 score of 0\.851, corresponding to only a modest improvement of 0\.012 over XGBoost\. The primary aim of this study, however, was not to optimize an operational tree species mapping approach, but to investigate which spectral bands and temporal periods contribute to model predictions and class separation\.

Combining recurrent models with Integrated Gradients enabled attribution for each individual spectral band for each available timestep without prior temporal aggregation\. In addition to conventional class\-wise attributions, using observations from an opposing class as the baseline allowed us to identify the features driving class separation, rather than only those contributing to the prediction of an individual class\. For pure and mixed Norway spruce, the contributions increased for B6, B7, B8, and B8A after approximately DOY 125, while B11 showed particularly strong attribution around DOY 100\. For beech and oak, B6 and B11 and the period between approximately DOY 110 and 150 were especially relevant for class separation\. The analysis further demonstrated that features contributing strongly to the prediction of an individual class are not necessarily identical to those driving the distinction between similar classes\.

The spectral temporal attribution patterns identified for pure and mixed Norway spruce broadly corresponded with known phenological processes, particularly the development of broadleaf vegetation within mixed stands\. The explainability analysis therefore provides a physically plausible interpretation of the model process, although it did not reveal substantially different spectral\-temporal relationships beyond those expected from existing phenological knowledge\. Still, explainable AI implications likeIGimprove the understanding of the discrimination process and can inform future features or data design\.

The increased number of annual observations offered by Sentinel\-2C motivates repeating the proposed workflow with denser time series in future work, as more observations improved classification performance in this study\. Additional temporal information may also reveal spectral\-temporal relationships that are currently difficult to resolve\. Furthermore, moving beyond categorical pure\-versus\-mixed labels toward estimating species proportions or subpixel forest composition could provide a more continuous representation of mixed stands and allow a more detailed investigation of the spectral response to increasing species admixture\. Overall, this study demonstrates that recurrent deep learning models combined with class\-comparative attribution provide a framework for investigating spectral\-temporal decision processes in forest classification, particularly where classes are spectrally and phenologically similar\.

## References

- \[1\]\(2024\)A Review: Tree Species Classification Based on Remote Sensing Data and Classic Deep Learning\-Based Methods\.Forests15\(5\),pp\. 852\(en\)\.External Links:ISSN 1999\-4907,[Link](https://www.mdpi.com/1999-4907/15/5/852),[Document](https://dx.doi.org/10.3390/f15050852)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[2\]J\. Stoffels, S\. Mader, J\. Hill, W\. Werner, and G\. Ontrup\(2012\)Satellite\-based stand\-wise forest cover type mapping using a spatially adaptive classification approach\.European Journal of Forest Research131\(4\),pp\. 1071–1089\(en\)\.External Links:ISSN 1612\-4677,[Link](https://doi.org/10.1007/s10342-011-0577-2),[Document](https://dx.doi.org/10.1007/s10342-011-0577-2)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[3\]J\. Tan, J\. Li, T\. Ma, X\. Yan, and Z\. Huo\(2025\)Leveraging Sentinel\-1/2 time series and deep learning for accurate forest tree species mapping\.Frontiers in Forests and Global Change8,pp\. 1599510\.External Links:ISSN 2624\-893X,[Link](https://www.frontiersin.org/articles/10.3389/ffgc.2025.1599510/full),[Document](https://dx.doi.org/10.3389/ffgc.2025.1599510)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[4\]S\. I\. Seneviratne, X\. Zhang, M\. Adnan, W\. Badi, C\. Dereczynski, A\. Di Luca, S\. Ghosh, I\. Iskandar, J\. Kossin, S\. Lewis, F\. Otto, I\. Pinto, M\. Satoh, S\. M\. Vicente\-Serrano, M\. Wehner, and B\. Zhou\(2021\)Weather and Climate Extreme Events in a Changing Climate\.InClimate Change 2021: The Physical Science Basis\. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change,V\. Masson\-Delmotte, P\. Zhai, A\. Pirani, S\. L\. Connors, C\. Péan, S\. Berger, N\. Caud, Y\. Chen, L\. Goldfarb, M\. I\. Gomis, M\. Huang, K\. Leitzell, E\. Lonnoy, J\. B\. R\. Matthews, T\. K\. Maycock, T\. Waterfield, O\. Yelekçi, R\. Yu, and B\. Zhou \(Eds\.\),pp\. 1513–1766\.External Links:[Document](https://dx.doi.org/10.1017/9781009157896.013)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[5\]J\. Bauhus, D\. I\. Forrester, and H\. Pretzsch\(2017\)Mixed\-Species Forests: The Development of a Forest Management Paradigm\.InMixed\-Species Forests,H\. Pretzsch, D\. I\. Forrester, and J\. Bauhus \(Eds\.\),pp\. 1–25\(en\)\.External Links:ISBN 9783662545515 9783662545539,[Link](http://link.springer.com/10.1007/978-3-662-54553-9_1),[Document](https://dx.doi.org/10.1007/978-3-662-54553-9%5F1)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1),[§1](https://arxiv.org/html/2610.19270#S1.p2.1)\.
- \[6\]S\. Yachi and M\. Loreau\(1999\)Biodiversity and ecosystem productivity in a fluctuating environment: The insurance hypothesis\.Proceedings of the National Academy of Sciences96\(4\),pp\. 1463–1468\(en\)\.External Links:ISSN 0027\-8424, 1091\-6490,[Link](https://pnas.org/doi/full/10.1073/pnas.96.4.1463),[Document](https://dx.doi.org/10.1073/pnas.96.4.1463)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[7\]European Commission\. Joint Research Centre\.\(2020\)Climate change impacts and adaptation in Europe: JRC PESETA IV final report\.\.Publications Office,LU\(eng\)\.External Links:[Link](https://data.europa.eu/doi/10.2760/171121)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[8\]A\. Iglseder, M\. Immitzer, A\. Dostálová, A\. Kasper, N\. Pfeifer, C\. Bauerhansl, S\. Schöttl, and M\. Hollaus\(2023\)The potential of combining satellite and airborne remote sensing data for habitat classification and monitoring in forest landscapes\.International Journal of Applied Earth Observation and Geoinformation117,pp\. 103131\(en\)\.External Links:ISSN 15698432,[Link](https://linkinghub.elsevier.com/retrieve/pii/S1569843222003193),[Document](https://dx.doi.org/10.1016/j.jag.2022.103131)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[9\]M\. Martone, P\. Rizzoli, C\. Wecklich, C\. González, J\. Bueso\-Bello, P\. Valdo, D\. Schulze, M\. Zink, G\. Krieger, and A\. Moreira\(2018\)The global forest/non\-forest map from TanDEM\-X interferometric SAR data\.Remote Sensing of Environment205,pp\. 352–373\(en\)\.External Links:ISSN 00344257,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0034425717305795),[Document](https://dx.doi.org/10.1016/j.rse.2017.12.002)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[10\]M\. Immitzer, F\. Vuolo, and C\. Atzberger\(2016\)First Experience with Sentinel\-2 Data for Crop and Tree Species Classifications in Central Europe\.Remote Sensing8\(3\),pp\. 166\(en\)\.External Links:ISSN 2072\-4292,[Link](https://www.mdpi.com/2072-4292/8/3/166),[Document](https://dx.doi.org/10.3390/rs8030166)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1),[§3\.1](https://arxiv.org/html/2610.19270#S3.SS1.p1.1)\.
- \[11\]D\. Phiri, M\. Simwanda, S\. Salekin, V\. Nyirenda, Y\. Murayama, and M\. Ranagalage\(2020\)Sentinel\-2 Data for Land Cover/Use Mapping: A Review\.Remote Sensing12\(14\),pp\. 2291\(en\)\.External Links:ISSN 2072\-4292,[Link](https://www.mdpi.com/2072-4292/12/14/2291),[Document](https://dx.doi.org/10.3390/rs12142291)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[12\]E\. Grabska, P\. Hostert, D\. Pflugmacher, and K\. Ostapowicz\(2019\)Forest Stand Species Mapping Using the Sentinel\-2 Time Series\.Remote Sensing11\(10\),pp\. 1197\(en\)\.External Links:ISSN 2072\-4292,[Link](https://www.mdpi.com/2072-4292/11/10/1197),[Document](https://dx.doi.org/10.3390/rs11101197)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[13\]L\. Blickensdörfer, K\. Oehmichen, D\. Pflugmacher, B\. Kleinschmit, and P\. Hostert\(2024\)National tree species mapping using Sentinel\-1/2 time series and German National Forest Inventory data\.Remote Sensing of Environment304,pp\. 114069\(en\)\.External Links:ISSN 00344257,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0034425724000804),[Document](https://dx.doi.org/10.1016/j.rse.2024.114069)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1),[§1](https://arxiv.org/html/2610.19270#S1.p2.1),[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p1.1),[§5\.3](https://arxiv.org/html/2610.19270#S5.SS3.p3.1)\.
- \[14\]M\. Wegler, P\. Kacic, F\. Thonfeld, S\. Holzwarth, N\. Jaggy, U\. Gessner, and C\. Kuenzer\(2025\)Tree species from space: a new product for Germany based on Sentinel\-1 and \-2 time series\.International Journal of Remote Sensing46\(16\),pp\. 6075–6108\(en\)\.External Links:ISSN 0143\-1161, 1366\-5901,[Link](https://www.tandfonline.com/doi/full/10.1080/01431161.2025.2530236),[Document](https://dx.doi.org/10.1080/01431161.2025.2530236)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1),[§1](https://arxiv.org/html/2610.19270#S1.p2.1),[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p1.1),[§3\.1](https://arxiv.org/html/2610.19270#S3.SS1.p1.1),[§5\.2](https://arxiv.org/html/2610.19270#S5.SS2.p2.1),[§5\.3](https://arxiv.org/html/2610.19270#S5.SS3.p3.1)\.
- \[15\]P\. Liu, C\. Ren, Z\. Wang, M\. Jia, W\. Yu, H\. Ren, and C\. Xia\(2024\)Evaluating the Potential of Sentinel\-2 Time Series Imagery and Machine Learning for Tree Species Classification in a Mountainous Forest\.Remote Sensing16\(2\),pp\. 293\(en\)\.External Links:ISSN 2072\-4292,[Link](https://www.mdpi.com/2072-4292/16/2/293),[Document](https://dx.doi.org/10.3390/rs16020293)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1),[§1](https://arxiv.org/html/2610.19270#S1.p2.1),[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p1.1)\.
- \[16\]A\. M\. Abdi and F\. Wang\(2026\)Mapping forest tree species and their uncertainty using Earth observation and National Forest Inventory data: towards operational monitoring in Sweden\.International Journal of Remote Sensing47\(7\),pp\. 2912–2943\(en\)\.External Links:ISSN 0143\-1161, 1366\-5901,[Link](https://www.tandfonline.com/doi/full/10.1080/01431161.2026.2625513),[Document](https://dx.doi.org/10.1080/01431161.2026.2625513)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1)\.
- \[17\]F\. Mouret, D\. Morin, M\. Planells, and C\. Vincent\-Barbaroux\(2025\)Tree Species Classification at the Pixel Level Using Deep Learning and Multispectral Time Series in an Imbalanced Context\.Remote Sensing17\(7\),pp\. 1190\(en\)\.External Links:ISSN 2072\-4292,[Link](https://www.mdpi.com/2072-4292/17/7/1190),[Document](https://dx.doi.org/10.3390/rs17071190)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1),[§1](https://arxiv.org/html/2610.19270#S1.p2.1)\.
- \[18\]M\. Freudenberg, S\. Schnell, M\. Henning, and P\. Magdon\(2025\)Consistent Multi\-Temporal National Scale Tree Species Maps for Germany 2018 – 2024\.External Links:[Link](https://www.ssrn.com/abstract=5705254),[Document](https://dx.doi.org/10.2139/ssrn.5705254)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p1.1),[§3\.1](https://arxiv.org/html/2610.19270#S3.SS1.p1.1),[§5\.3](https://arxiv.org/html/2610.19270#S5.SS3.p3.1)\.
- \[19\]H\. Abdullah, A\. K\. Skidmore, R\. Darvishzadeh, and M\. Heurich\(2019\)Timing of red\-edge and shortwave infrared reflectance critical for early stress detection induced by bark beetle \(Ips typographus, L\.\) attack\.International Journal of Applied Earth Observation and Geoinformation82,pp\. 101900\.External Links:ISSN 15698432,[Document](https://dx.doi.org/10.1016/j.jag.2019.101900)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p2.1)\.
- \[20\]M\. Freudenberg, S\. Schnell, and P\. Magdon\(2025\)A Sentinel\-2 machine learning dataset for tree species classification in Germany\.Earth System Science Data17\(2\),pp\. 351–367\(en\)\.External Links:ISSN 1866\-3516,[Link](https://essd.copernicus.org/articles/17/351/2025/),[Document](https://dx.doi.org/10.5194/essd-17-351-2025)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p2.1),[§2\.1](https://arxiv.org/html/2610.19270#S2.SS1.p1.1),[§3\.1](https://arxiv.org/html/2610.19270#S3.SS1.p2.1),[§3\.3](https://arxiv.org/html/2610.19270#S3.SS3.p1.1),[§5\.4](https://arxiv.org/html/2610.19270#S5.SS4.p1.1)\.
- \[21\]E\. u\. H\. Bundesministerium für Landwirtschaft\(2022\)Bundeswaldinventur 2022\.External Links:[Link](https://www.thuenen.de/de/themenfelder/waelder/die-bundeswaldinventur/ergebnisse-der-bundeswaldinventur-2022)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p2.1)\.
- \[22\]M\. Immitzer, M\. Neuwirth, S\. Böck, H\. Brenner, F\. Vuolo, and C\. Atzberger\(2019\)Optimal Input Features for Tree Species Classification in Central Europe Based on Multi\-Temporal Sentinel\-2 Data\.Remote Sensing11\(22\),pp\. 2599\(en\)\.External Links:ISSN 2072\-4292,[Link](https://www.mdpi.com/2072-4292/11/22/2599),[Document](https://dx.doi.org/10.3390/rs11222599)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p2.1),[§5\.2](https://arxiv.org/html/2610.19270#S5.SS2.p2.1)\.
- \[23\]A\.M\.C\. Malkanthi, R\.M\.K\.T\. Rathnayaka, B\.T\.G\.S\. Kumara, and D\.R\. Welikanna\(2025\)Forest Cover Classification with Explainable AI: A Comparative Analysis\.In2025 9th SLAAI International Conference on Artificial Intelligence \(SLAAI\-ICAI\),Colombo, Sri Lanka,pp\. 1–6\.External Links:ISBN 9798331561239,[Link](https://ieeexplore.ieee.org/document/11318453/),[Document](https://dx.doi.org/10.1109/SLAAI-ICAI68534.2025.11318453)Cited by:[§1](https://arxiv.org/html/2610.19270#S1.p2.1)\.
- \[24\]Thünen\-Institut\(2024\)Bundeswaldinventur 2022 – Ergebnisdatenbank: Vorrat nach Land und Baumartengruppe\.External Links:[Link](https://bwi.info/inhalt1.3.aspx?Text=3.03%20Baumartengruppe&prRolle=public&prInv=BWI2022&prKapitel=3.03)Cited by:[§2\.1](https://arxiv.org/html/2610.19270#S2.SS1.p1.1)\.
- \[25\]S\. J\. Russell and P\. Norvig\(2020\)Artificial intelligence: a modern approach\.Fourth Edition edition,Pearson Series in Artificial Intelligence,Pearson,Hoboken, NJ\.External Links:ISBN 978\-0\-13\-461099\-3Cited by:[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p1.1)\.
- \[26\]M\. Kirsch, J\. Wernicke, P\. S\. Datta, and C\. Preisach\(2025\)Early Detection of Forest Disturbances in Homogeneous Stands – Deep Learning\-Based Monitoring ofPicea Abies\.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing18,pp\. 22300–22316\.External Links:ISSN 1939\-1404, 2151\-1535,[Link](https://ieeexplore.ieee.org/document/11141690/),[Document](https://dx.doi.org/10.1109/JSTARS.2025.3603094)Cited by:[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p1.1),[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p3.1)\.
- \[27\]Y\. Kong, Q\. Huang, C\. Wang, J\. Chen, J\. Chen, and D\. He\(2018\)Long Short\-Term Memory Neural Networks for Online Disturbance Detection in Satellite Image Time Series\.Remote Sensing10\(3\),pp\. 452\.External Links:ISSN 2072\-4292,[Document](https://dx.doi.org/10.3390/rs10030452)Cited by:[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p1.1)\.
- \[28\]S\. Hochreiter and J\. Schmidhuber\(1997\)Long Short\-Term Memory\.Neural Computation9\(8\),pp\. 1735–1780\.External Links:ISSN 0899\-7667, 1530\-888X,[Document](https://dx.doi.org/10.1162/neco.1997.9.8.1735)Cited by:[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p1.1)\.
- \[29\]K\. Cho, B\. Van Merrienboer, C\. Gulcehre, D\. Bahdanau, F\. Bougares, H\. Schwenk, and Y\. Bengio\(2014\)Learning Phrase Representations using RNN Encoder–Decoder for Statistical Machine Translation\.InProceedings of the 2014 Conference on Empirical Methods in Natural Language Processing \(EMNLP\),Doha, Qatar,pp\. 1724–1734\(en\)\.External Links:[Link](http://aclweb.org/anthology/D14-1179),[Document](https://dx.doi.org/10.3115/v1/D14-1179)Cited by:[§A\.1](https://arxiv.org/html/2610.19270#A1.SS1.p4.1),[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p1.1)\.
- \[30\]J\. Lee, Y\. Lee, J\. Kim, A\. R\. Kosiorek, S\. Choi, and Y\. W\. Teh\(2018\)Set Transformer: A Framework for Attention\-based Permutation\-Invariant Neural Networks\.arXiv\.External Links:[Link](https://arxiv.org/abs/1810.00825),[Document](https://dx.doi.org/10.48550/ARXIV.1810.00825)Cited by:[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p3.1)\.
- \[31\]A\. Vaswani, N\. Shazeer, N\. Parmar, J\. Uszkoreit, L\. Jones, A\. N\. Gomez, L\. Kaiser, and I\. Polosukhin\(2017\)Attention Is All You Need\.arXiv\.External Links:[Link](https://arxiv.org/abs/1706.03762),[Document](https://dx.doi.org/10.48550/ARXIV.1706.03762)Cited by:[§2\.2](https://arxiv.org/html/2610.19270#S2.SS2.p4.1)\.
- \[32\]W\. Samek, T\. Wiegand, and K\. Müller\(2017\)Explainable Artificial Intelligence: Understanding, Visualizing and Interpreting Deep Learning Models\.arXiv\.External Links:[Link](https://arxiv.org/abs/1708.08296),[Document](https://dx.doi.org/10.48550/ARXIV.1708.08296)Cited by:[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p1.1)\.
- \[33\]V\. Hassija, V\. Chamola, A\. Mahapatra, A\. Singal, D\. Goel, K\. Huang, S\. Scardapane, I\. Spinelli, M\. Mahmud, and A\. Hussain\(2024\)Interpreting Black\-Box Models: A Review on Explainable Artificial Intelligence\.Cognitive Computation16\(1\),pp\. 45–74\(en\)\.External Links:ISSN 1866\-9956, 1866\-9964,[Link](https://link.springer.com/10.1007/s12559-023-10179-8),[Document](https://dx.doi.org/10.1007/s12559-023-10179-8)Cited by:[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p1.1)\.
- \[34\]M\. Ancona, E\. Ceolini, C\. Öztireli, and M\. Gross\(2019\)Gradient\-Based Attribution Methods\.InExplainable AI: Interpreting, Explaining and Visualizing Deep Learning,W\. Samek, G\. Montavon, A\. Vedaldi, L\. K\. Hansen, and K\. Müller \(Eds\.\),pp\. 169–191\.External Links:ISBN 978\-3\-030\-28954\-6,[Link](https://doi.org/10.1007/978-3-030-28954-6_9),[Document](https://dx.doi.org/10.1007/978-3-030-28954-6%5F9)Cited by:[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p1.1)\.
- \[35\]M\. Sundararajan, A\. Taly, and Q\. Yan\(2017\)Axiomatic Attribution for Deep Networks\.arXiv\.External Links:[Link](https://arxiv.org/abs/1703.01365),[Document](https://dx.doi.org/10.48550/ARXIV.1703.01365)Cited by:[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p1.1),[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p2.1),[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p6.1)\.
- \[36\]A\. Jha, J\. K\. Aicher, M\. R\. Gazzara, D\. Singh, and Y\. Barash\(2020\)Enhanced Integrated Gradients: improving interpretability of deep learning models using splicing codes as a case study\.Genome Biology21\(1\),pp\. 149\(en\)\.External Links:ISSN 1474\-760X,[Link](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02055-7),[Document](https://dx.doi.org/10.1186/s13059-020-02055-7)Cited by:[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p1.1)\.
- \[37\]S\. Sanyal and X\. Ren\(2021\)Discretized Integrated Gradients for Explaining Language Models\.InProceedings of the 2021 Conference on Empirical Methods in Natural Language Processing,Online and Punta Cana, Dominican Republic,pp\. 10285–10299\(en\)\.External Links:[Link](https://aclanthology.org/2021.emnlp-main.805),[Document](https://dx.doi.org/10.18653/v1/2021.emnlp-main.805)Cited by:[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p1.1)\.
- \[38\]L\. Tronchin, E\. Cordelli, L\. R\. Celsi, D\. Maccagnola, M\. Natale, P\. Soda, and R\. Sicilia\(2024\)Translating Image XAI to Multivariate Time Series\.IEEE Access12,pp\. 27484–27500\.External Links:ISSN 2169\-3536,[Link](https://ieeexplore.ieee.org/document/10439172/),[Document](https://dx.doi.org/10.1109/ACCESS.2024.3366994)Cited by:[§2\.3](https://arxiv.org/html/2610.19270#S2.SS3.p1.1)\.
- \[39\]B\. Schieber, R\. Janík, and Z\. Snopková\(2009\)Phenology of four broad\-leaved forest trees in a submountain beech forest\.Journal of Forest Science55\(1\),pp\. 15–22\.External Links:ISSN 12124834, 1805935X,[Link](http://jfs.agriculturejournals.cz/doi/10.17221/51/2008-JFS.html),[Document](https://dx.doi.org/10.17221/51/2008-JFS)Cited by:[§5\.2](https://arxiv.org/html/2610.19270#S5.SS2.p1.1)\.
- \[40\]K\. E\. Lewińska, O\. Hagolle, F\. Gascon, P\. Griffiths, U\. Leser, and P\. Hostert\(2026\)Benefits of the Sentinel\-2 mission triplet constellation – one year on\.Remote Sensing Applications: Society and Environment44,pp\. 102278\(en\)\.External Links:ISSN 23529385,[Link](https://linkinghub.elsevier.com/retrieve/pii/S2352938526004118),[Document](https://dx.doi.org/10.1016/j.rsase.2026.102278)Cited by:[§5\.3](https://arxiv.org/html/2610.19270#S5.SS3.p2.1)\.
- \[41\]R\. Verdecchia, L\. Cruz, J\. Sallou, M\. Lin, J\. Wickenden, and E\. Hotellier\(2022\)Data\-Centric Green AI An Exploratory Empirical Study\.In2022 International Conference on ICT for Sustainability \(ICT4S\),Plovdiv, Bulgaria,pp\. 35–45\.External Links:ISBN 9781665482868,[Link](https://ieeexplore.ieee.org/document/9830097/),[Document](https://dx.doi.org/10.1109/ICT4S55073.2022.00015)Cited by:[§5\.3](https://arxiv.org/html/2610.19270#S5.SS3.p3.1)\.

## Appendix AAppendix

### A\.1Equations for LSTM and GRU

ft=σ⁡\(wf​xt\+uf​ht−1\+bf\)f\_\{t\}=\\sigma\\left\(w\_\{f\}x\_\{t\}\+u\_\{f\}h\_\{t\-1\}\+b\_\{f\}\\right\)\(4\)it=σ⁡\(wi​xt\+ui​ht−1\+bi\)i\_\{t\}=\\sigma\\left\(w\_\{i\}x\_\{t\}\+u\_\{i\}h\_\{t\-1\}\+b\_\{i\}\\right\)\(5\)ot=σ⁡\(wo​xt\+uo​ht−1\+bo\)o\_\{t\}=\\sigma\\left\(w\_\{o\}x\_\{t\}\+u\_\{o\}h\_\{t\-1\}\+b\_\{o\}\\right\)\(6\)C~t=t​a​n​h​\(wC~​xt\+uC~​ht−1\+bC~\)\\tilde\{C\}\_\{t\}=tanh\\left\(w\_\{\\tilde\{C\}\}x\_\{t\}\+u\_\{\\tilde\{C\}\}h\_\{t\-1\}\+b\_\{\\tilde\{C\}\}\\right\)\(7\)Ct=ft⊙Ct−1\+it⊙C~tC\_\{t\}=f\_\{t\}\\odot C\_\{t\-1\}\+i\_\{t\}\\odot\\tilde\{C\}\_\{t\}\(8\)ht=ot⊙t​a​n​h​\(Ct\)h\_\{t\}=o\_\{t\}\\odot tanh\\left\(C\_\{t\}\\right\)\(9\)
Where:

- •The*forget gate*is denotedftf\_\{t\}, the*input gate*iti\_\{t\}, the*cell candidate gate*C~t\\tilde\{C\}\_\{t\}and the*output gate*oto\_\{t\}\.
- •The input weights are represented asww, the recurrent weights for the Hidden State asuuand the biases are denoted asbb\.
- •hhis the hidden state,CCis the cell state andxxis the input data at time steptt\.
- •⊙\\odotdenotes the point\-wise multiplication of two vectors\.
- •σ\\sigmais the sigmoid activation function,*tanh*is the tanh activation function\.

For the GRU unit the gates are defined by the following functions\[[29](https://arxiv.org/html/2610.19270#bib.bib8)\]:

rt=σ⁡\(wr​xt\+br\+wr​ht−1\+br\)r\_\{t\}=\\sigma\\left\(w\_\{r\}x\_\{t\}\+b\_\{r\}\+w\_\{r\}h\_\{t\-1\}\+b\_\{r\}\\right\)\(10\)zt=σ⁡\(wz​xt\+bz\+wz​ht−1\+bz\)z\_\{t\}=\\sigma\\left\(w\_\{z\}x\_\{t\}\+b\_\{z\}\+w\_\{z\}h\_\{t\-1\}\+b\_\{z\}\\right\)\(11\)gt=tanh⁡\(wg​xt\+bg\+rt⊙\(wg​ht−1\+bg\)\)g\_\{t\}=\\tanh\\left\(w\_\{g\}x\_\{t\}\+b\_\{g\}\+r\_\{t\}\\odot\\left\(w\_\{g\}h\_\{t\-1\}\+b\_\{g\}\\right\)\\right\)\(12\)ht=\(1−zt\)⊙gt\+zt⊙ht−1h\_\{t\}=\\left\(1\-z\_\{t\}\\right\)\\odot g\_\{t\}\+z\_\{t\}\\odot h\_\{t\-1\}\(13\)
Where:

- •The*reset gate*is denoted byrtr\_\{t\}, the*update gate*byztz\_\{t\}, and the*candidate hidden state*bygtg\_\{t\}\.

The sigmoid activation function is calculated as:

σ⁡\(x\)=11\+e−x\\sigma\(x\)=\\frac\{1\}\{1\+e^\{\-x\}\}\(14\)
The tanh activation function is calculated as:

t​a​n​h​\(x\)=ex−e−xex\+e−xtanh\(x\)=\\frac\{e^\{x\}\-e^\{\-x\}\}\{e^\{x\}\+e^\{\-x\}\}\(15\)

### A\.2Additional Tables

Table 5:Comparison of model configurationsTable 6:Hyperparameter search spaces for Random Forest andXGBoostmodelsModelHyperparameterTested valuesRandom Forestn\_estimators100, 250, 500, 750max\_depth3, 6, 10, 15min\_samples\_split2, 5, 10min\_samples\_leaf1, 2, 5max\_featuressqrt,log2, 0\.5, 0\.7bootstrapTruemax\_samples0\.5, 0\.7, 0\.9,Noneclass\_weightNone,balancedXGBoostmax\_depth3, 6, 10, 15learning\_rate0\.01, 0\.1, 0\.2, 0\.3, 0\.4subsample0\.5, 0\.6, 0\.7, 0\.8, 0\.9colsample\_bytree0\.5, 0\.6, 0\.7, 0\.8, 0\.9colsample\_bylevel0\.5, 0\.6, 0\.7, 0\.8, 0\.9n\_estimators100, 250, 500, 750reg\_alpha0\.1, 0\.001, 0\.00001reg\_lambda0\.1, 0\.001, 0\.00001RNN\-ModelsRNNLSTM, GRUAttentionMHA, NoneBi\-directionalTrueBottleneck dimension32, 64, 128Hidden DimensionBottleneck dimension×4\\times 4dropout0\.2