The study area (approx. 35° 58′–37° 2′ E and 111° 45′–112° 32′ N) is located in the Shanxi Province, China. It includes the Huodong National Coal Mining Area and Taiyue Mountain National Forest Park, the largest forest reservation in Shanxi (Figure 1). The Taiyue Mountain and surrounding forests form an important nature reserve with a temperate continental climate, an average annual temperature of 9.2 °C, and 564 mm of precipitation. The study area covers a total of 7,715.17 km², with elevations from 534 to 2,564 meters, and forests occupy 3,861.14 square kilometers, approximately 50% of the region, which belongs to the temperate northern forest zone. The main dominant tree species include: Larix principis-rupprechtii (LP), Pinus tabuliformis (PT), Pinus bungeana (PB), Platycladus orientalis (PO), Quercus wutaishanica (QW), Betula spp. (BA), Populus spp. (PS).

Based on the GEE platform, we used Sentinel-1 Ground Range Detected (GRD) backscatter products and Sentinel-2 MultiSpectral Instrument (MSI) Level-2A products from 2022 to 2023, including cloud removal, monthly median value extraction, and resampling to a spatial resolution of 10 meters, visit https://developers.google.com/earth-engine/datasets/catalog/sentinel for more details. We calculated commonly used vegetation indices and additional indices recognized in prior researches as valuable for vegetation remote sensing with Sentinel (Frampton et al., 2013; Li et al., 2014; Ngo et al., 2023; Schulz et al., 2024). Combined indices with the original bands, 33 variables include: B2-B8, B8A, B11-12, VV, VH, NDVI, GNDVI, LSWI, EVI, NDVIre1-3, NDre1-2, Clre, PSRI, MSAVI, MSRre, MTCI, CCCI, S2REP, RVI, VVVHR, NDIVV. Full names and formulas can be found in Supplementary Table S1. Temporal profiles of these indices for different tree species are provided in Supplementary Figures S1, S2, the plots highlight species-specific spectral, phenological, and structural differences, such as distinct temporal patterns between deciduous and evergreen trees, and backscatter contrasts between coniferous and broadleaf trees (Supplementary Figure S3).

The European Space Agency (ESA) WorldCover 10 m v200 product was used in this study to extract forest areas mask.

As reference, we used the Forest Inventory data provided by Taiyue Mountain forest administration. The data was collected around 2020 and 2021 through field surveys and high-resolution images, recording forest resource attributes such as dominant and secondary species, stand age and area in polygons. We selected seven main tree species in the study area as target classes, categorizing other species as “Others” (mostly broadleaf) for the final tree species mapping (Figure 1). Reference points was randomly generated within each polygon, imposing constraints of over 20 meters between points and 30 meters from plot edges, and the number of points was based on plot area, resulting in a reference dataset. To address the simplification of mixed forest scenarios and accuracy overestimation resulting from the use of only pure-species samples, we also selected mixed-species plots to create pseudo-labeling samples, combining them with pure samples for DL model training and evaluation. The specific methods are detailed in section 3.1. We employed stratified sampling to select 60% of the samples from each species as an independent testing dataset to prevent data leakage (Table 1). The remaining 40% of samples were used for training and validation with five-fold cross-validation, each fold involved using one subset for validation. This process ensured validation on entirely different samples, aiming to thoroughly evaluate model performance, ensure generalization ability, and adjust model hyperparameters.

Additionally, forest mask was used to generate unlabeled pretraining forest samples by creating grid points at 100 m intervals (10 rows/columns for images), resulting in over 380,000 unlabeled samples.

Forest inventory data consists of polygonal plots documenting the dominant and secondary tree species. Previous studies commonly used samples from pure units for model training and accuracy assessment (Immitzer et al., 2019; Yang et al., 2024), which may reduce dataset representativeness, lead to lower model performance in identifying mixed-species areas and overestimate mapping accuracy (Blickensdörfer et al., 2024). It’s easier to distinguish between evergreen and deciduous species due to their greater phenological feature differences. To ensure credibility for the pseudo samples added to the dataset, only mixed units of evergreen and deciduous species were used. Specifically, pure units of target tree species were selected firstly to generate pure samples, then evergreen-deciduous mixed units among the seven target tree species were selected to generate unlabeled samples. The process for generating pseudo-labels for unlabeled samples based on deep learning classification model are as Figure 2. First, the pretrained model was fine-tuned using labeled samples from pure units of two species (one evergreen and one deciduous), to obtain a binary classification model. Next, the unlabeled samples of mixed units were input into the binary classification model, which outputs predicted scores ranging from 0 to 1 for each sample. The class with the highest score was assigned as the pseudo-label, and samples with highest scores above 0.9 were added to the final dataset.

To further analyze the forest ecosystem in the study area, we calculated biodiversity indices based on the final tree species classification map: Species Richness, the Shannon-Wiener Index, and Simpson’s Diversity Index (Peng et al., 2021). Species Richness reflects diversity by simply counting the number of species in a sample. The Shannon-Wiener Index takes into account both species richness and evenness, based on the relative abundance of each species in the sample. Simpson’s Diversity Index emphasizes dominant species by considering the square of relative abundance, highlighting species evenness. The formulas are as follows:

Where R is Species Richness, H is Shannon-Wiener Index, D is Simpson’ s Diversity Index, S is total number of species in sample, and p_i is the proportion of individuals of species i relative to the total number of individuals. we calculated each index at a resolution of 100 m.

We generated a confusion matrix from validation results and calculated six commonly used metrics: Overall Accuracy (OA), Kappa, User’s Accuracy (UA), Producer’s Accuracy (PA), F1 Score and macro-F1. UA (Precision) and PA (Recall) focus on class-specific accuracy, while the F1 Score combines both for a comprehensive measure. OA, Kappa and macro-F1 evaluate overall performance.

Where k is the number of classes, TP is the number of true positives, FP is the number of false positives, FN is the number of false negatives,N_correct is the number of corrected classified samples, N_total is the number of all samples, P_o is the overall accuracy, P_e is the proportion of agreement expected by chance.

Considering the potential bias caused by the spatial distribution of samples, we adopted an inverse sampling-intensity weighted method (De Bruin et al., 2022). This method corrects estimation bias by assigning more weight to sparsely sampled areas and less weight to densely sampled areas, based on the sample distribution density. Using the Kernel Density Estimation (KDE) method from the scikit-learn package in Python (Pedregosa et al., 2011), we estimated the sample point density and applied inverse sampling weights. Additional accuracy metrics, including weighted-OA, weighted-kappa, weighted-F1 and weighted-macro-F1, were calculated to provide a more comprehensive evaluation of accuracy.

Five-fold cross-validation was employed to evaluate the model performance and generalization ability based on 2 years data, with accuracy variations and F1-scores for each species shown in Figure 4. The model showed stable performance across different folds, with an average Kappa of 0.78, OA of 0.82, and macro-F1 of 0.80. The F1-scores for each tree species were also consistent, with average values as follows: LP 0.83, PT 0.88, PB 0.76, PO 0.83, QW 0.76, PA 0.84, PS 0.76, OT 0.75. The lowest standard deviation was for PT 0.003, while the LP, PB, PO, QW, and BA were around 0.013, highest standard deviation was for OT 0.037. These results indicate that the model is stable and reliable for tree species classification, and the final accuracy to be confirmed by the testing samples.

We first tested the model using only data from June to July 2022, and the results achieved an OA of only 0.49, Kappa of 0.38 and macro-F1 of 0.38. We then gradually extended the time series by adding 1 month of data at both the start and end in each iteration, analyzing the impact of varying time series length (Figure 5). The results showed that increasing the input length gradually improves the performance of the model. When the series was extended to include data from April to September, which is commonly used for vegetation and crop remote sensing studies, the classification accuracy significantly improved, with an OA of 0.68, Kappa of 0.62, and macro-F1 of 0.64. This time range covers the greenness rise and fall period for most tree species, as well as the leaf-on and off period for deciduous species. As the temporal depth increased to cover the late leaf-fall period of deciduous trees, feature differences between evergreen and deciduous species grew, resulting in a gradual improvement in each classification accuracy metrics. The OA approached 0.8, with a macro-F1 of 0.78, using monthly data for all of 2022, demonstrating effective species differentiation.

Additionally, we incorporated data from the following year (2023) into the time series, adding one month at a time (13–24 months represent January to December 2023). This approach aimed to determine whether extending the time series beyond an entire growth cycle would enhance accuracy of the model and capture more potential time dependencies and hidden features across annual data. The results showed that extending the length slightly improved classification accuracy, but the gains were much smaller than those seen when the temporal data did not cover a complete growth cycle. Moreover, simply extending the time series does not guarantee improved accuracy. For instance, when the data was extended by an additional 3 months (i.e., 1–15 months) beyond the full year of 2022, here was little improvement in accuracy, with the macro-F1 fluctuating around 0.77 and no increase in OA or Kappa. This indicating that the new data may lack significant additional information. As data from the second year (2023) growing season was incorporated, all metrics continued to grow gradually. Extending the temporal length to 24-months, covering 2 years, all accuracy metrics reached their highest level: OA 0.847, Kappa 0.815, and macro-F1 0.836, Both the accuracy metrics and the F1 score growth curve approached a stable state (Figure 5). We also tested the model using only Sentinel-2 24-months data for classification, and the results showed slightly lower accuracy compared to the combination of two data sources, with an OA of 0.819, Kappa of 0.782, and macro-F1 of 0.805 (confusion matrix in Supplementary Table S2).

The Inverse sampling-intensity weighted method was used to account for potential estimation bias from the spatial distribution of samples. Additional weighted accuracy metrics were calculated, with the confusion matrix shown in Table 2 (Estimated sampling intensity and sample weights distribution in Supplementary Figure S4). The results show a slight decrease in F1 scores for several species, with PO and BA decreasing from 0.85 and 0.86 to 0.81, respectively. There were also adjustments in the overall accuracy metrics, with OA decreasing from 0.847 to 0.834 and macro-F1 dropping from 0.836 to 0.813. These adjustments further improved the reliability of the results, and the overall performance remains satisfactory.

By using 24 months data, the model demonstrated its ability to produce high-quality tree species classification. We applied the trained model to process time series data of all forest pixels, producing a tree species distribution map for the entire region (Figure 6). The map showed good spatial consistency with the forest inventory reference data, demonstrating the reliability of the classification outcomes. The area statistics for each tree species are as follows: Pinus tabuliformis covers the largest area at 2,368.54 km², followed by Quercus wutaishanica at 760.14 km². Larix principis-rupprechtii, Pinus bungeana and Betula spp. have similar coverage, with 138.64 km², 145.14 km², and 136.29 km², respectively. Populus spp. occupies a small area of 17.36 km², while Others (mainly broadleaf) cover 123.53 km².

Three biodiversity indices were calculated using the tree species distribution map and displayed the results (Figure 7). The Species Richness ranges from 1 to 8, indicating the number of tree species within a unit. The Shannon-Wiener Index ranges from 0 to 1.949, with higher values indicating a more complex and even tree species composition. The Simpson’s Diversity Index ranges from 0 to 0.851, with values closer to 1 indicating that most individuals are concentrated in a few species, reflecting lower evenness. These results illustrate the biodiversity of the study area forest, offering insights into the health and resilience of ecosystem.

The global separability index (SI-global) combines the feature variable separability of various tree species, as shown in Figure 8A, where the intensity of colors indicates differentiation of a feature variable among tree species at a specific time point, reflecting its overall capability to distinguish between species. The results showed that the vegetation indices (e.g., NDVI, GNDVI, NDre1, CIre, etc.) calculated from NIR and red-edge bands demonstrate higher separability, compared to the original band values from satellite sensors. The specific separability indices of all pairs are shown in Supplementary Figure S5. We found that the distinction between deciduous and evergreen species exhibited significantly separability, and extracted the separability index map (Figures 8B,C) for four pairs of tree species (one pair of evergreen species (PT-PO), one pair of deciduous species (LP-BA), and two pairs of mixed deciduous and evergreen species (PT-BA and LP-PT).) The SI heatmaps for same-type species pairs showed low separability, indicating minimal differences in shallow features. In contrast, different species pairs exhibited significant separability, particularly from January to April, October to December, and July to September, with higher index values during these periods. Most deciduous trees are in leafless stages during these periods, leading to significant differences from evergreens in spectral reflectance. These results also highlighted the challenge of distinguishing species with similar phenology.

PCA was used for dimensionality reduction on the dataset with 33 variables. The results (Figure 9A) showed that the majority of the variance in tree species samples is concentrated in PC1 (37.91%) and PC2 (19.83%). We further analyzed the composition of the two principal components by displaying the distribution of variables (Figure 9B) and the top 10 important variables for each component (Figures 9C,D). Vegetation indices have high contributions in PC1, while optical bands dominate PC2, and SAR data shows minimal contribution in the PCA analysis. In PC1, the red edge indices (NDre1, NDre2, MSReN, MSRre, NDVIre1, Clre) and the optical vegetation indices calculated using near-infrared bands (NDVI, mSAVI, GNDVI, LSWI) contribute significantly, these variables primarily reflect the physiological indicators such as chlorophyll. PC2 is primarily influenced by the original optical satellite bands, reflecting the surface forest canopy’s spectral information, such as leaf color and brightness. To further investigate the importance of each variable at different time points, PCA was performed on monthly data (2023) to assess the contributions to the first two PCs, ranking changes are shown in Supplementary Figure S6, and monthly variable loadings in Supplementary Figure S7. While contributions fluctuated over time, the trend remained consistent with previous analysis. Vegetation indices made the largest contribution to PC1, consistently showing the highest contributions across different time. In contrast, PC2 was primarily influenced by the Sentinel-2 bands.

The pretrain and fine-tune approach was employed to enhance the classification model ability to identify tree species. We trained a model without pretraining, and compared the changes in OA and macro-F1 during the training process with pretrained model, as well as the final F1 scores for each tree species (Figure 13). The results demonstrated that the model, which was pretrained on forest pixels, significantly improved the accuracy for classification, with noticeable enhancements in the identification accuracy across all species. Additionally, the accuracy curve over training epochs shows that the pretrained model reached the performance ceiling of the non-pretrained model around the 20th epoch, further boosting the accuracy ceiling. The cross-validation results further confirm the enhanced consistency of the pretrained model (Supplementary Figure S9). Existing studies have shown that pretraining models can effectively enhance generalization ability and performance (Jing and Tian, 2021; Misra and Van Der Maaten, 2020). But creating a pretraining task that aligns well with the target task is challenging and requires further research, such as transfer learning strategies and general models (Miller et al., 2024).

Extending time series data generally improved classification accuracy for tree species in most cases (Figure 5), highlighting the advantage of using time series images with more varying information (Foerster et al., 2012). Using only June and July data resulted in an OA of just 0.49. Despite the period being the peak growing season for most plants, the shorter series lacked phenological variation, leaving the model to rely almost entirely on spectral reflectance values. Trees with similar phenological traits tend to exhibit highly similar spectral signatures, but key phenological changes are crucial for species identification, and longer time series capturing these shifts provide more information for better differentiation (Liu et al., 2023). Extending the data beyond a full year can still slightly improve accuracy of identifying tree species, but the gains were less significant than adding within-year data, and some additional data may not contribute valuable information. Visualizing the model soft outputs showed that adding an extra year of data boosts prediction confidence, allowing better identification of subtle differences between tree species and optimizing predictions (Figure 12). As time series length increases, classification accuracy for individual tree species and overall performance stabilized at a satisfactory level, and we believe that further extension may not significantly enhance true classification accuracy.

The physiological differences between deciduous and evergreen trees result in distinct spectral reflectance patterns, leading to a more pronounced separability between these two groups. The results indicated that during the leaf fall period, distinct types tree species show significant separability (Figure 8). In contrast, similar growth processes among the same type species lead to insignificant difference in feature values across time series, making accurate identification difficult using only original bands and indices. Many studies emphasize the significance of NIR and Red-edge indices for vegetation research (Frampton et al., 2013; Li et al., 2014), we also found that original bands are less effective for distinguishing between species compared to index variables, even when comparing two distinct types of trees. Vegetation indices effectively highlight differences between tree species by sensitively reflecting physiological indicators like chlorophyll and water content. PCA results also confirm the importance of vegetation index variables (Figure 9), with a high contribution from vegetation indices in the most significant principal component (PC1). PC1 primarily captures seasonal subtle variations in tree physiological characteristics, while PC2, influenced mainly by original bands, reflects differences in leaf color, consistent with previous studies (Schulz et al., 2024). There is variability in the importance of different variables to the principal components across different time periods, due to differences in sensitivity to phenological characteristics at various stages. However, the primary contributors to the two principal components remain the vegetation indices and spectral bands, respectively. We also calculated Pearson’s correlation coefficients among the 33 feature variables to assess inter-variable relationships (Supplementary Figure S10). The results revealed higher correlations and possible redundancy among optical bands with similar wavelengths and their derived vegetation indices, as well as among SAR-derived variables. In contrast, the correlations between variables from the two different sources were generally low. This suggests that adding SAR data to optical inputs may provide complementary information for distinguishing tree species, as the combined data achieved higher tree species classification accuracy compared to using only optical features, with macro-F1 scores of 0.836 and 0.805, respectively. This aligns with previous studies integrating multiple data sources, suggesting that additional information provided by SAR, such as vegetation structural characteristics, may contribute to improved classification performance, which is also reflected in the distinct backscatter profiles across species (Supplementary Figure S2).

In addition to producing high-accuracy tree species classification maps using deep learning models on SITS, visualizing the model can help us understand its internal processes and bolster its reliability. Multi-head attention weight matrices revealed the model behavior in capturing information from various aspects of the input time series, including key phenological changes and periods of features similarity (section 4.5.1). The model employed nonlinear transformations and multi-head self-attention mechanisms to continuously extract and integrate features, enhancing class-relevant information while diminishing irrelevant features. This process evolves the input features into a more abstract and separable form, illustrating the significant advantage of deep learning in processing time series data and capturing information (Vaswani et al., 2017; Xu et al., 2020).

Projecting hidden features from each model layer and the original features, clearly illustrated how low-level inputs are transformed into high-level representations that effectively distinguish tree species (section 4.5.2). In the shallow layers of model, input retains low-level features closely tied to their original forms, leading to mixed distributions of different tree species. As the model progresses to deeper layers, it captures features of increasing complexity most relevant to the categories, transforming hidden features to project different tree species into a more distinct feature space for precise separation (Goodfellow et al., 2016). We also found that the samples of deciduous species cluster closely after projection, as do the evergreen species, which aligns with our feature analysis, indicating smaller differences and reduced separability among similar tree species. In the final output layer, we observed significant inter-class separation among samples, indicating that deeper phenological features are captured by the model. Similarities lead to close proximity in the projected space, while differences create distinct clusters.

The tree species distribution map and forest biodiversity results revealed a certain correlation between tree species distribution and elevation (Figure 6). Pinus tabuliformis dominates the forests of study area, primarily located in the low-elevation regions of eastern Taiyue Mountain and near the Huodong coal mining area, which also exhibited lower levels of biodiversity. According to information gathered from forestry departments, this trend may be attributed to large-scale plantings in the 1980s and reclamation plantings after mining activities, leading to a lack of species diversity. In contrast, Taiyue Mountain is predominantly covered by primary forests, which exhibited a richer species composition and higher levels of biodiversity. Many studies have found that primary forests have higher species diversity and biomass than planted and second-growth forests (Cavanaugh et al., 2014; Shirima et al., 2015). The second largest tree species in the study area, Quercus wutaishanica, is mainly in the southern part of Taiyue Mountain. Pinus bungeana is distributed along the western edge of Taiyue, forming a continuous strip, while Platycladus orientalis is located in the northwest, near the Pinus bungeana. Larix principis-rupprechtii and Quercus wutaishanica are concentrated in the high-altitude central region. The distribution and biodiversity data of tree species provided crucial information and effective support for forest ecosystem management and ecological research. For mining areas, where forest harvesting and replanting is a recurring process, it is crucial to implement more comprehensive and diverse planting strategies to protect and enhance both forest biodiversity and carbon sequestration capacity.