Assessing Landsat-8 and Sentinel-2 spectral-temporal features for mapping tree species of northern plantation forests in Heilongjiang Province, China

Mengyu Wng, Yi Zheng, Chengqun Hung, Rn Meng, Yong Png, Wen Ji,Jie Zhou, Zehu Hung, Linhun Fng, Feng Zho,＊

a Key Laboratory of Geographical Process Analysis & Simulation of Hubei Province/College of Urban and Environmental Sciences, Central China Normal University,Wuhan, 430079, China

b Department of Geographical Sciences, University of Maryland, College Park, MD, 20742, USA

c Macro Agriculture Research Institute, Interdisciplinary Sciences Research Institute, College of Resources and Environment, Huazhong Agricultural University, Wuhan,430070, China

d Institute of Forest Resource Information Techniques, Chinese Academy of Forestry, Beijing, 100091, China

e Key Laboratory of Forestry Remote Sensing and Information System, National Forestry and Grassland Administration, Beijing, 100091, China

Keywords:Tree species mapping Plantation forests Red-edge features Temporal frequency of data acquisition Fusion of Landsat-8 and Sentinel-2

ABSTRACT

1. Introduction

The planted forests provide important ecosystem services and represent valuable global resources, such as carbon sequestration and timber products (Wingfield et al., 2015; Chen et al., 2019; FAO, 2021). Plantation forest is a particular type of planted forests where planted trees are managed intensively and usually composed of one or two even-aged species with regular spacing (FAO, 2021). As home to the world's largest plantation forests, China has contributed to about one fourth of net global leaf area increase in the past two decades through land-use management activities such as afforestation, especially in northern China(Chen et al.,2019;FAO,2021).Common timber species in China's northern plantation forests include larch, Mongolian pine, Korean pine and Mongolian oak, etc., which are also keystone species in Eurasian temperate-boreal ecozones (Hytteborn et al., 2005). These northern plantation forests provide valuable timber resources and support regional economic development, but their ecosystem services have been threatened by climate change and increasing disturbances (Hanewinkel et al.,2013; Wingfield et al., 2015; Harris et al., 2021). Accurate mapping of these tree species is essential for effective management of northern plantation forests in China(Zhang et al.,2000;Meng et al.,2022)and is also valuable for characterizing ecosystem services and forest-climate feedbacks in this region (Bonan et al., 1992; Duveiller et al., 2021).Traditional forest surveys could produce detailed and accurate maps of tree species distribution, but they are time-consuming and labor-intensive, and subjecting to artificial errors that could not be traced.

Remote sensing has emerged as one of the most efficient tools for mapping tree species at regional (Lim et al., 2019; Kollert et al., 2021)and landscape scales (Shen and Cao, 2017; Fang et al., 2020; Lindberg et al., 2021). Because of its ability to penetrate the canopy and characterize tree structure, airborne Lidar data can successfully differentiate tree species and produce highly accurate species maps (Dalponte et al.,2012;Blomley et al.,2017;M¨ayr¨a et al.,2021).However,airborne Lidar data are often expensive and have limited spatial coverage, making it challenging for large-scale applications(Heinzel and Koch,2011;Li et al.,2013). Spaceborne Lidar data, such as data from the ICESat or GEDI missions, are more cost-effective, but they do not provide enough point cloud densities for distinguishing tree species at fine spatial scales(Queinnec et al.,2021;Silva et al.,2021;Leite et al.,2022).Radar data,especially L-band or P-band SAR, also provide helpful structural information to improve tree species mapping accuracies (Wong and Fung,2014;Dost'alov'a et al.,2021;Zhao et al.,2022).But the inherent speckle noises in SAR data pose huge challenges for extracting structural features from the data and it is not feasible to use SAR data alone for accurate mapping of tree species (Bjerreskov et al., 2021). As a result, for large-scale tree species mapping applications, optical remote sensing features (i.e., spectral and temporal features) are more commonly used than structural information from Lidar or SAR data.

For remote sensing-based tree species mapping,spectral information is key for distinguishing different species. Hyperspectral images were found effective in mapping tree species in many ecosystems (Dalponte et al., 2013; Shang and Chisholm, 2014; Meng et al., 2018a, 2018b;Wietecha et al.,2019;Modzelewska et al.,2020;Zhang et al.,2020;Yang et al.,2020),but the high costs and limited availability of high-resolution hyperspectral data hinders its application in fine-scale tree species monitoring that is relevant to forest management. Compared with natural forests, plantation forests are generally composed of even-aged single species and have simpler forest structure and species composition(Wingfield et al.,2015;Yu et al.,2020). Consequently, the spectral separability among plantation tree species is generally higher than that of the natural forests, making it possible to distinguish different species with multi-spectral measurements.

Recently,temporal features of multi-spectral imagery received much attention for tree species mapping. Approaches using time series multispectral images, mostly spaceborne, were found useful for mapping tree species(Immitzer et al.,2019;Hemmerling et al.,2021;Chen et al.,2021).There are also some recent developments in tree species mapping using multi-spectral or multi-temporal images from airborne or UAV platforms(Cao et al.,2016;Shen and Cao,2017;Xu et al.,2020;Egli and H¨opke, 2020; Johansen et al., 2020; Belcore et al., 2021; Grybas and Congalton,2021). Despite of their high spatial resolutions,images from airborne or UAV platforms are limited in its spatial coverage and it is impractical to acquire repeated airborne observations for large areas.On the contrary, satellite-borne multi-spectral images often have coarser spatial resolutions than airborne images, but readily available globally and repeatedly.Some images from recent commercial satellites also have very high spatial resolutions (i.e., sub-meter) and can be used for mapping tree species at fine spatial scales (Pu and Landry, 2012; Immitzer et al.,2012;Waser et al.,2014;Li et al.,2019;Ferreira et al.,2019;Meng et al.,2017).But the high costs of these high resolution satellite images limit their applications in tree species mapping.

Many people are turning to free satellite images with medium resolution, such as Landsat images, for solution. Although the spatial and spectral resolutions of Landsat image are not high,the temporal sequence of Landsat images will likely favor the discrimination of different tree species,as each species may show distinct phenology across the growing season and time series Landsat observations will likely capture this spectral-temporal separability (Pasquarella et al., 2018; Hemmerling et al.,2021).Therefore,the potential of these time series Landsat images in classifying tree species in plantation forests worth further investigations(Rocchini et al.,2016).Among all the methods for utilizing spectral-temporal information from time series Landsat data, harmonic analysis is one of the most popular ones, given its ability to capture the intra-annual variations among different vegetation types and conditions(Nguyen et al., 2020; Roy and Yan,2020). This method can reconstruct time series feature of vegetation phenology from limited number and irregular time of observations, mostly via NDVI time series, meanwhile eliminating and replacing cloud-contaminated values (Menenti et al.,1993; Zhou et al., 2016). NDVI is most often used in harmonic analysis because it is a good indicator for vegetation phenology (Zhou et al.,2016). Harmonic regression coefficients from time series Landsat data were found to outperform median composite components in predicting forest canopy cover in conterminous USA(Derwin et al.,2020).Tasseled Cap harmonic parameters derived from all available Landsat imagery(1985–2015) were found useful for mapping general forest types (e.g.,hardwood swamp, softwood swamp, etc.), and the spectral-temporal approach outperformed the single and multiple image approaches (Pasquarella et al., 2018). Despite the substantial contributions of these studies, a knowledge gap still existed about the effects of temporal frequency of data acquisitions on the spectral-temporal approach for mapping tree species in plantation forests.

In addition to the time series Landsat images, the availability of Sentinel-2 (S2) imagery since year 2015, with its fine spatial, spectral,and temporal resolutions, provides more possibilities for tree species mapping. S2 sensors (including S2A and S2B) were designed to obtain multispectral images at finer spatial and temporal resolutions than Landsat (Astola et al., 2019; Clark, 2020). The newly added red-edge bands of S2 sensors are sensitive to variations of chlorophyll contents in the leaves and have been shown to improve the monitoring capability of plant nutrition(Westergaard-Nielsen et al.,2021),health status(Zhen et al.,2021)and vegetation mapping(Immitzer et al.,2016,2019).The increased temporal resolution and spectral bands of S2 imagery is promising for improving tree species and forest attributes mapping over Landsat-8(L8), yet there were some discrepancies among findings from different studies.Korhonen et al.(2017)compared the performance of S2 and L8 in estimating boreal forest canopy cover and leaf area index,and found that S2 performed marginally better than L8. Chrysafis et al.(2017)assessed S2 and L8 for mapping forest stock volume and obtained almost the same accuracies from these images. In a tree species classification test in a Mediterranean natural forest, Clark (2020) found that predictor variables derived from S2 (excluding red-edge related variables) were significantly better than comparable L8 variables.Despite the valuable insights provided by these studies,a comprehensive evaluation is still needed to examine the key spectral domains of L8 and S2 in mapping tree species in northern plantation forests,and to understand whether and to what extent the S2 red-edge bands can improve forest species mapping in these regions.

Harmonizing L8 and S2 images can further increase the temporal frequency of data acquisition and was found helpful in improving the accuracies of some vegetation monitoring and land cover classification studies than using single sensor alone (Carrasco et al., 2019; Griffiths et al., 2019). Although it is a general assumption that with increasing number of images,time-series based classification is prone to producing spectral redundancies and provide little additional information about plant condition, the correlations between S2 and L8 bands may vary in time-depending on the respective image acquisition frequency and the phenology of the depicted trees-and the mentioned assumption may not necessarily hold true for all tree species(Pe～na and Brenning,2015;Pe～na et al.,2017;Pu,2021).Some studies evaluated the performance of L8 and S2 fusion in mapping crop types and characterizing forest attributes(Liu et al., 2018, 2020; Htitiou et al., 2021). L8 and S2 together provided adequate observations to extract phenological characteristics for mapping sugarcane plantation (Wang et al., 2020). Fusion of L8 and S2 enabled improved mapping of cropping intensity in major crop producing areas in China(Liu et al.,2020).However,few studies have evaluated the effects of increased temporal frequency of data acquisition (from 16-day to 5-day and shorter) for tree species mapping in plantation forests.A better understanding of the separate and combined capabilities of these two sensors helps to improve tree species mapping in plantation forests over large spatial-temporal domains and thus is highly valuable.

In this study,we compared the utility of L8 and S2 in mapping several keystone plantation tree species(e.g.,larch,Mongolian pine,and Korean pine) in northern plantation forests of China. Specific aims of the study are:

1) Identify the key spectral domains of L8 and S2 for separating several keystone species in northern plantation forests across different phenological stages(greenness rise,greenness fall,leaf on and leaf off stages)as well as for the whole time-series,with a particular interest in the contributions of S2 red-edge features.

2) Evaluate whether, and to what extent, increasing the temporal frequency of data acquisitions(i.e.,single image,multi-date,time series of L8, S2, and L8 and S2 fusion) can improve tree species mapping accuracies in northern plantation forests.

2. Materials and methods

2.1. Study area and tree species reference data

Our study area locates in Heilongjiang, the most northeastern province of China sharing a border with Russia(Fig.1).This region is home to the largest forest ecosystem in China and produces one third of the country's timber product each year, including some of the keystone species of the temperate-boreal ecozones of the northern hemisphere(i.e.,larch,Mongolian Pine,Korean Pine).Forestry in this area has been an important contributor to regional economy since 1950s. The Mengjiagang forest is one of the largest and well-managed plantation forests in Heilongjiang Province with a total management area of about 155 km2(dark green areas in Fig.1),of which 135 km2are plantation forest.We chose Mengjiagang plantation forest as our study area because of its representativeness of northern plantation forests in China.It includes the common keystone species in the temperate-boreal ecozones in northern China (e.g., pure stands of larch, Mongolian Pine, Korean Pine, Mongolian Oak,and the mixtures of these species).As one of the most important plantation forests in northern China,the Mengjiagang plantation also has accurate reference data so that we can fully evaluate the L8 and S2 spectral and temporal features of each species and forest type.

Fig.1. Study area of the Mengjiagang forest plantation.Dark green areas show the management area of the plantation and colored dots show the general distribution of the major timber species and mixed forest types.The left-bottom inset show the location of the study area in China.The base image is a Sentinel-2 image acquired on May 22,2020,displayed with the following band combination:red,Band 4;green,Bands 3;blue,Bands 2.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

In this study, we conducted our classification experiments based on the above tree species and forest types with abundant samples - all undisturbed forest areas during 2016–2020 as our samples,which resulted in 131,246 samples at 30-m resolution(Table 1). Beginning and ending dates for key phenological stages of these species were identified from literature and confirmed with local forestry experts (Luan et al., 1992;Guo et al., 2011; Cheng and Sun, 2013; Zhang, 2013) (Fig. 2). We acquired the wall-to-wall sub-compartment map produced by the Mengjiagang plantation administration during the forest resource investigation in year 2016 (Fig. S2). The sub-compartment map was a shapefile, generated by manually drawing of polygons according to the high spatial resolution orthophoto aerial imaging in year 2016,and then modified in year 2017 during field visit.For each polygon,there include attributes such as species,stand age,area,etc.

2.2. Data preparation and test setups for tree species classification

We have designed three sets of classification experiments - singledate, multi-date and spectral-temporal analysis (Fig. 3). For the singleand multi-date analysis, we selected 7 pairs of clear L8 and S2 observations obtained within two days apart between 2016 and 2020(Table 2).For the spectral-temporal analysis,we used all cloud-free L8 and S2 images for the study area during 2016–2020 to ensure temporal consistency with our training and validation samples (see more details of these images in supplementary material Fig. S1). All the L8 and S2 imagery was terrain-corrected and atmospherically corrected to surface reflectance.All S2 images were resampled to 30-m resolution using the bilinear interpolation, and the spectrum was resampled to be consistent with L8 following Zhang et al. (2018). Cloud and shadow of both L8 and S2 images were masked using the cloud removal algorithm implemented in Google Earth Engine (GEE). Then the preprocessed L8 and S2 images were downloaded from GEE.We first examined the key spectral domains for separating these tree species across different phenological stages(single-date experiment,Fig.3).The contributions of S2 red edge bands were examined by comparing the performances of S2 with and without red edge bands during key phenological stages.

To enhance the discrimination power of spectral features, multiple spectral bands are usually synthesized into vegetation indices (VIs). In this study, we extracted the surface reflectance and common VIs corresponding to the samples for L8 images and the spatially and spectrally resampled S2 images (supplementary material Table S1 and S2). For single-date classifications, there were 32 features for L8 or S2 without red-edge and 48 features for S2 with red-edge. For multi-date classification,we combined spectral features from two images(with the highest accuracies) from each of the phenological stages (Fig. 3). Then in the spectral-temporal experiments, we also compared the accuracies of S2 with and without red edge temporal features. We only changed the red edge related features and all other features remained constant.

To extract spectral-temporal features from time-series datasets, we fitted and calculated the harmonic coefficients (amplitude, phase andconstant) from each sample. As the harmonic analysis algorithm is sensitive to noise, we manually selected images with good quality to construct the time series and extract temporal features. Three sets of experiments with different temporal frequency of data acquisition were used as inputs: (1) sparse L8 time series (14 images), (2) dense S2 time series(34 images),and(3)super dense L8 and S2 fusion time series(48 images).We then fitted all valid data at pixel scale with the normal leastsquares regression model(Eq.(1)).

Table 1 The four tree species and three mixed forest type classes mapped in this study.Class sample locations were extracted from the local forest resource survey map.

We used higher frequencies harmonic model to capture more temporal variations from the time series. Each band or index corresponded five harmonic parameters. To reduce computation burden and avoid over-fitting, we limited the spectral-temporal features to the harmonic coefficients extracted from the original bands and the most important VIs selected by VSURF (Variable Selection Using Random Forest) from the single-date experiments (see section 2.3 for details).

To explore the effects of different temporal frequency of data acquisition for tree species classification, we prepared three datasets, with increasing temporal frequency of data acquisitions, as inputs for classification experiments and feature selection: (1) single-date L8 and S2 images,(2)multi-date L8 and S2 images,and(3)L8 time series,S2 time series and the fusion of L8 and S2 time series images(Fig.3).The benefits of S2 red edge bands were examined for all three datasets. For each dataset,we extracted L8 and S2 sample features from GEE,referencing to the forest resource survey map.To coincide with the reference data,both L8 and S2 images between January 1,2016 and December 31,2020 were used in this study. To minimize the effects of disturbances on our classification experiments, we have applied the LandTrendr algorithm implemented in GEE (Kennedy et al., 2018) to remove the disturbed forest areas between 2016 and 2020.

2.3. Tree species classification and feature selection

All classification experiments were performed using the Random Forest (RF) algorithm (Breiman, 1999) implemented in the R package randomForest. We set ntree (the number of random trees) to 500 and tuned mtry(the number of input variables used at each node)from 2 to k,where k is the number of features. The samples were randomly divided into a training set and a test set, equaling to 70% and 30% of the total samples,respectively.

Feature selection is vital for reducing the dimension of the highdimensional datasets and preventing model overfitting (Grabska et al.,2020), and it is also helpful for us to understand the most important spectral bands and indices for tree species classification.In this study,we used VSURF procedure(Genuer et al.,2016)for feature selection,which is an R algorithm based on the variable importance metric of random forest classification. VSURF returns two feature subsets: an important subset with some redundancy and a smaller one trying to eliminate redundancy (Genuer et al., 2016; Grabska et al., 2020). We used the second subset from VSURF as our final classification inputs (see parameters we used in supplementary material Table S4).

2.4. Accuracy assessments

Based on the test samples shown in Table 1, we evaluated the classification accuracies, including the overall accuracy (OA, Eq. (2)), producer's accuracy(PA,Eq.(3)),user's accuracy(UA,Eq.(4))and macro-F1 score (Eq. (5)), for all the experiments. We also produced classification maps to examine the spatial consistency with the forest resource survey map,and calculated areas of each tree species and forest types.

Fig. 2. The phenological calendar for major timber tree species in the study area.

Fig. 3. Workflow for this study. VSURF: Variable Selection Using Random Forest.

where TP was true positive, representing the number of correct predictions for positive samples; FN was false negative, indicating the number of incorrect predictions for negative samples; TN was true negative, indicating the number of correct predictions for negativesamples; FP was false positive, representing the number of incorrect predictions for positive samples; x was the class index.

Table 2 Date of the Landsat-8 and Sentinel-2 images used for the single-date experiment.The image dates obtained from the two sensors were within two days apart to minimize the influence of phenological changes.

3. Results

3.1. Spectral signatures of major tree species

The spectral curves of the four tree species showed the highest separability in the NIR region during greenness rise and leaf on stages,in the SWIR during greenness fall, and in the visible region during leaf off(Fig.4).During greenness rise and leaf on stages,the Mongolian oak had the highest NIR reflectance (for both B8 and B8A), followed by larch,Korean pine and Mongolian pine. S2 red edge bands (especially B7)showed similar separability as the NIR band, while other spectrum showed little separability among the species.For the greenness fall stage,Mongolian oak also showed the highest SWIR reflectance, followed by larch, Mongolian pine and Korean pine. During the leaf-off stage, the Mongolian oak still had the highest reflectance,especially in the visible bands, followed by larch, Mongolian pine and Korean pine. For all the phenological stages, the red edge bands all showed similar or less separability than the NIR bands. For mixed forest types, the highest separability between mixed broadleaf forest and Mongolian oak were in the NIR and red edge bands during greenness rise,SWIR in the greenness fall,and visible bands in leaf off. The highest separability between mixed broadleaf forest and larch was consistent with the above. Mixed conifer and Mongolian pine as well as Korean pine can be most easily identified in the NIR and visible bands during the greenness fall and leaf off stages,respectively.

Spectral signatures from S2 showed higher separability than that of the L8, even for very similar types such as Mongolian oak and mixed broadleaf, indicating improved spectral performances of S2 sensor over L8. The additional red edge bands and narrow band NIR of S2 consistently showed similar or less separability than the broadband NIR band for all the phenological stages.

Fig. 4. Snapshots of average Landsat-8 and Sentinel-2 spectra for the four timber tree species and three mixed forest types during different phenological stages.

The spectral-temporal features (e.g., NDVI time series) of different species also showed distinct patterns(Fig. 5). Evergreen needleleaf species, Mongolian pine and Korean pine, had the flattest temporal curve and the NDVI values of Korean pine were consistently higher than that of Mongolian pine across the annual cycle.For deciduous species,broadleaf Mongolian oak showed larger variations of NDVI values than narrowleaf larch,with both higher peak NDVI and lower bottom NDVI.

3.2. Comparing L8 and S2 for tree species mapping

S2 showed better performances than L8 in all the classification experiments, even without the red edge bands (Fig. 6). The accuracy differences between L8 and S2 were highest in May(3.3%and 4.6%for OA and macro-F1, respectively) and lowest in October (0.3% and 0.5% for OA and macro-F1, respectively). Regarding PA and UA, classification accuracies based on S2 images were better than L8 in almost all cases,especially for PA during greenness rise,leaf on and leaf off stages(Fig.7).PAs were lowest for mixed conifer forests and mixed conifer and broadleaf and using S2 images can improve the PA for mixed conifer and broadleaf for about 20%. S2 red edge bands did not show significant improvements in species mapping accuracies in the classification experiments.

Fig. 5. Observed (dots) and fitted (lines) temporal NDVI values, as inputs for time-series harmonic analysis, for the four tree species and three mixed forest types. Harmonic coefficients (amplitude, phase and constant) were further extracted from the fitted time-series NDVI curve as spectral-temporal features for tree species classifications.

NDTI(the ratio of SWIR1 minus SWIR2 to SWIR1 plus SWIR2)and the Tasseled Cap coefficients(i.e.,TCW and TCG)were selected as important spectral features for distinguishing four timber tree species and three mixed forest types throughout key phenological stages(Table 3).During the leaf on stage, NIR and SWIR1 were important spectral features, and during the leaf off stage,NDTI and the visible bands associated features(i.e., NormG, ARVI, RGRI and NormGR) played important roles. Even though red edge features were selected by VSURF as important features in single-date, multi-date and S2 spectral temporal experiments, especially during the greenness rise and fall stages, no additional accuracy gains were achieved by adding the red edge features.

For separating tree species within the deciduous forests, NDTI was one of the most important spectral features throughout all phenological stages(Table 4).In transition(greenness rise and greenness fall)and leaf on stages,the indices associated with NIR band(i.e.,LSWI,DVI,SARVI,GNDVI and RVI) were also important. In leaf off stage, the indices associated with visible bands outperformed other indices. In distinguishing evergreen classes, the red edge related features showed outstanding separability in transition(greenness rise and greenness fall)and leaf off stages,especially the RE1 related indices.The most important spectral features for separating tree species within evergreen classes during the leaf on stage was TCW.

3.3. Effects of temporal frequency of data acquisition for tree species mapping

Fig. 6. Overall accuracy and macro-F1 score for species classifications. S–T: spectral-temporal; Fusion: fusion of Landsat-8 and Sentinel-2 without red edge. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

In China's northern plantation forests, with relatively homogenous forest stands and simple species composition, increasing the temporal frequency of data acquisition can increase the accuracies for tree species mapping for up to 3.2%(from 90.1%using single-date imagery to 93.3%using S2 spectral-temporal features,Fig.6).For single date analysis,the highest accuracies were achieved using S2 images from the greenness rise and fall images (90.1% and 89.9% when using S2 images obtained in May and October, respectively), likely due to the maximized spectral differences of evergreen needleleaf (i.e., Mongolian pine and Korean pine), deciduous needleleaf (i.e., larch), and deciduous broadleaf (i.e.,Mongolian oak) species during sprouting and senescence. Multi-date classification showed improved classification results compared with single-date classification, especially when combining complementary leaf on image and greenness fall or greenness rise images (92.9% and 92.6% when using S2 images obtained in June plus October and June plus May, respectively). Use of spectral-temporal images provided slightly better classification accuracies than using multi-date images,which is the case for both S2 and L8. Accuracies from S2 time series(93.3%) were marginally better than that of the L8 time series (93.0%)and the best multi-date accuracy,and adding L8 to S2 time series did not improve classification accuracies (93.2%), indicating that for China's northern plantation forests, the effects of increasing the temporal frequency of data acquisition could saturate quickly after using only two images from key phenological stages.

According to the feature selection results of three time series groups(L8 time series,S2(without red edge features)time series and the fusion of L8 and S2 time series), the spectral-temporal features of the VIs associated with the red bands (i.e., ARVI, NDSVI, RGRI, etc.) showed important contributions for distinguishing major tree species in China's northern plantation forests. Time-series normalized difference VIs related to red edge bands and narrow NIR (i.e., NDVIre1n, NDVIre2n,NDVIre3n) as well as RE1-related indices (i.e., CIre, MSRren) were also found as useful features.

3.4. Forest species map and area statistics

The predicted tree species map showed good spatial consistency with the forest resource survey map(Fig.8).Mixed broadleaf and Mongolian oak took over most of the Northeastern corner of the plantation while the other species spread across the study area in a patchy manner. Mixed broadleaf forest was the largest forest type (Fig. 9). Larch, as the dominant timber species,occupied more than 20 km2of the plantation.The second largest timber species was the Mongolian pine, followed by Korean pine, Mongolian oak, and mixed conifer. The area of mixed broadleaf and conifer was less than 5 km2.

4. Discussions

Results from our study show that with near-equivalent bands and reduced spatial resolution, S2-based tree species mapping outperforms that of L8 (0.4%–3.4% and 0.2%–4.4% higher for OA and macro-F1,respectively) in tree species mapping for northern plantation forests in China,especially during the greenness rise and greenness fall stages.This might be explained by the fact that S2 is down-sampled from 10 to 30 m resolution and includes more spatial information for separating some of the species, such as Mongolian pine and mixed conifer (Figs. 4 and 7).This finding is consistent with some previous studies about the utilities of S2 for forest applications.Topaloˇglu et al.(2016)and Astola et al.(2019)pointed out that even when S2 is down-sampled to the same spatial resolution as L8,the spatial information carried by S2 pixels still provides additional information for forest classification. Besides the rich spatial information that S2 can provide,other potential reasons for the superior performance of S2 might be its narrower visible bands and broader NIR band,which are worth of exploration in future studies.

S2-based tree species mapping outperforms that of L8 across phenological stages in our study, with or without the red edge bands (Fig. 6).Previous studies show that although red edge bands can be helpful for forest attributes estimation and forest stress monitoring,they had limited or marginal capability in improving the accuracy metrics(Delegido et al.,2011; Korhonen et al., 2017; Zarco-Tejada et al., 2018; Bhattarai et al.,2020). Results from our feature selection analysis coincide with these studies:the red edge bands were selected as important features when the leaves are changing the most (e.g., greenness rise and greenness fall)(Table 3),but adding red edge bands do not result in significant accuracy gains. It is possible that for tree species classification in northern plantation forests,red edge bands do not provide additional information than other bands (e.g., visible, NIR and SWIR), who are also sensitive to changes in chlorophyll and leaf areas (Houborg et al., 2009). Thus,adding or removing red-edge bands does not significantly affect the classification outcomes in the northern plantation forests.

Fig. 7. Spider charts representing the producer's accuracy and user's accuracy for seven tree species in key phenological stages. LG: larch; PS: Mongolian pine; PK:Korean pine; QM: Mongolian oak; MC: mixed conifer; MB: mixed broadleaf; MCB: mixed conifer and broadleaf.

Our results demonstrate that the use of spectral-temporal features provided improved performances in tree species classification than using multi-date and individual images, which is the case for both S2 and L8.This likely can be explained by the fact that for different phenological stages, NIR, SWIR and the visible bands provided high separability for these species (Fig. 4), and the spectral temporal features can combine these spectral differences and maximize the separability, and thus outperform features from single or multiple images. Although we can also reconstruct time-series spectral features from any single year and map annual tree species, given that we achieved the highest accuracies using spectral-temporal features extracted from 5-year time series data,using images from a single year will likely reduce the mapping accuracy.This is because with much fewer observations,the harmonic analysis willnot be able to do as a good job in fitting the phenological trends of each species(Zhou et al.,2016;Pasquarella et al.,2018).

Table 3 Important spectral features for discriminating four tree species and three mixed forest types, selected by VSURF. The green shaded features are red edge bands related. Please see the full name, equation and reference for these features in the supplementary material.

Although some previous tree species classification research adopted the hierarchical classification method and achieved good classification results,we find that for separating tree species in our study area,a direct classification approach with the spectral-temporal features works the best (see results of our spectral-temporal approach and the hierarchical classification approach in supplementary material Tables S5 and S6,respectively). In a hierarchical classification, forest areas are often classified into evergreen and deciduous forest first,and then classify specific species within the evergreen and deciduous area masks (Chen et al.,2020;Illarionova et al.,2021).This approach is especially helpful when higher level class separability is much higher than lower class separability, and by dividing the multi-class classification task into several levels,it has more control over uncertainties associated with each level of classification (Sun and Lim, 2001). The downsides of the hierarchical classification method are error propagations among layers and complex processing procedures(Silla and Freitas,2010).In our study,we did not use the hierarchical classification method as these tree species have enough separability for a direct classification and the added steps and processing of the hierarchical classification did not result in improved classification accuracies.

For mapping plantation tree species in our study area, the model using spectral-temporal features from time series S2 achieved the highest classification accuracy and adding L8 images to S2 time series did not further improve S2 classification results (Fig. 6). These findings are consistent with results from a previous study on mapping natural forest types using Landsat spectral-temporal features(Pasquarella et al.,2018),except that for plantation forests,the contribution of temporal frequency of data acquisition can saturate more quickly than natural forests. The use of multi-date images can achieve nearly as good accuracies as the spectral-temporal feature-based classification, which are important information for making decisions regarding what data and methods to use for mapping tree species in this region or regions with similar climate and species composition.For example,if a contemporary tree species map is needed for similar studies areas,time series S2 and our selected spectral indices will be most helpful;or if time series species maps dating back to years without the S2 data are necessary,then the Landsat time series or even multi-date classification from Landsat could probably provide comparable classification accuracies.

Our results also showed that S2 time series provided similar accuracies as fusion of L8 and S2 time series(Fig.6).One possible explanation might be that the increased temporal frequency of data acquisition does not necessarily mean added information: affected by clouds and snow,the clear images of L8 and S2 we used to extract harmonic parameters all clustered in several months of year(Fig.S1).The combination of L8 and S2 might have limited contribution in improving the fitting of harmonic analysis and the extraction of the subtle temporal variations. More intelligent temporal data mining algorithms, such as LSTM (Hochreiter and Schmidhuber,1997)or NTM(Graves et al.,2014)from the recurrent neural network family, might be promising future directions for further improving tree species mapping accuracies using spectral-temporal features(Zhong et al.,2019).

The lowest classification accuracies (especially producer's accuracy)of the final classification results are from mixed conifer and broadleaf forest and mixed conifer forest(Fig.7 and Table S6).Some mixed conifer pixels were mistakenly classified as larch or Mongolian pine. According to the forest resource survey map,larch or Mongolian pine were often the dominant tree species in the mixed conifer forest type, making it challenging to separate mixed conifer and these pure stands. In addition,many mixed conifer and broadleaf pixels were mistakenly classified as the mixed broadleaf or larch.This is because mixed conifer and broadleaf forest usually include larch and other broadleaf species.Therefore,these classes are highly similar in both spectral and temporal features,leading to the over-estimation of larch and mixed broadleaf and underestimation of mixed conifer and mixed conifer and broadleaf (Fig. 9). Even in transitional stages,the separability between the mixed forests with larch are low. For future studies, images with higher spectral or spatial resolutions might be helpful for increasing the separability of these species and further improving mapping accuracies.

In this study, we have identified important spectral features for separating all tree species(Table 3)and species within the evergreen and deciduous classes (Table 4), consistent with the spectrum identified as having the highest spectral separability in Fig.4.These spectral featureswill be helpful for forest managers or researchers who would like to build their own tree species classification models from L8 or S2 in ecosystems with similar species compositions. In addition, we can apply the best performing model to a larger study area to support regional forest management,such as the Huanan county where Mengjiagang plantation forest locates in(shown in supplementary material Fig.S3).Although we include most of the key timber tree species in northern China, such as larch,Mongolian pine,Korean pine and Mongolian oak,other typical tree species in northern China, such as poplar (Populus girinensis), Chinese pine (Pinus tabulaeformis), and locust (Robinia pseudoacacia) were not included, and a more comprehensive evaluation is suggested for future studies.Since our study area has little forest change(less than 1%)during 2016–2020, we only focused on stable forests, and our mapping results can represent species distribution during the whole study interval. By excluding the changed areas for the study interval,we also minimized the uncertainties induced by the time difference between validation sample points (acquired in 2016) and the image acquisition (2016–2020). For areas with frequent forest disturbances, however, it will be helpful to extract pre-disturbance spectral-temporal features and identify the most commonly disturbed tree species by agents (i.e., fires, insects, etc.), to recognize priority species for future forest management (Rogers et al.,2015;Wingfield et al.,2015).

Table 4 Important spectral features for discriminating species within the deciduous and evergreen classes: (a) Larch and Mongolian oak; (b)Mongolian pine and Korean pine, selected by VSURF. The green shaded features are features based on the red edge bands.Please see the full name, equation and reference for these features in the supplementary material.

Our study area includes some of the keystone species in northern China and the temperate-boreal ecozone,such as larch,Mongolian pine and Mongolian oak. These species are widely used in plantation forests worldwide, providing important ecosystem services and pillaring regional economic development (Wingfield et al.,2015). Recent studies show that plantation forests are major carbon sinks around the world(Lei et al., 2019; Tong et al., 2020). Large-scale mapping of the species distribution and compositions are thus critical for quantitative evaluations of the carbon stored in these plantation trees.Accurate tree species maps are also critical for tracking species composition(Stanke et al.,2021)and scheduling harvesting activities(Ceccherini et al.,2020)to maximize the carbon sequestration potential of these plantation forests, in the face of climate change.

5. Conclusions

In this study,we comprehensively evaluated the contributions of key spectral domains of L8 and S2 and different temporal frequency of data acquisitions to mapping several keystone tree species in China's northern plantation forests.Results indicate that:1)S2 outperformed L8 images in all the tree species classification experiments, especially when using images from the greenness rise and greenness fall stages. With almost equivalent bands, S2 still significantly outperformed L8 classification accuracies, indicating that S2's superior performances over L8 in mapping key tree species in China's northern plantation forests were likely due to the improved spatial resolution, rather than the added red-edge bands; 2) In all the classification experiments, NDTI (the ratio of SWIR1 minus SWIR2 to SWIR1 plus SWIR2) and Tasseled Cap coefficients related features were among the most important features for distinguishing key tree species in China's northern plantation forests,and spectral-temporal features extracted from time-series, red band related VIs (i.e., ARVI, NDSVI, RGRI, etc.) were also useful for the time seriesbased classifications; 3) Use of spectral-temporal images provided slightly better classification accuracies than using multi-date images,which is the case for both S2 and L8. Accuracies from S2 time series(93.3%) were marginally better than that of the L8 time series (93.0%)and the multi-date accuracy(92.9%),and adding L8 to S2 time series did not improve classification accuracies(93.2%),indicating that for China's northern plantation forests, the effects of increasing the temporal frequency of data acquisition could saturate quickly after using only two images from key phenological stages. Including key tree species compositions in the temperate-boreal ecozone,this study thus has important implications for monitoring similar forest species composition over large scale and contributing to regional sustainable development.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Funding

This work was supported by National Natural Science Foundation of China (Grant No. 41901382), Open Fund of State Key Laboratory of Remote Sensing Science (Grant No. OFSLRSS201917), and the HZAU research startup fund(No. 11041810340;No.11041810341).

Fig. 8. Predicted tree species map in Mengjiagang plantation forest,with insets a–d showing the comparisons between the forest resource survey map(a–d(1))and the predicted map (a–d (2)) across the study area. This map was produced based on the best classification model results (i.e., S2 spectral-temporal classification).

Fig. 9. Areas of the four tree species and three mixed forest types in the study area. Reference areas are calculated from the forest resource survey map, while the predicted areas are from the best results of this study (i.e., S2 spectral-temporal classification).

Availability of data and materials

The data and materials generated for this study is available from the corresponding author on reasonable request.

Authors' contributions

Mengyu Wang and Yi Zheng performed the experiments,analyzed the data, prepared figures and tables, and prepared the first draft of the manuscript.Chengquan Huang,Ran Meng,Jie Zhou and Linchuan Fang edited and validated the manuscript with critical comments and reviewed the manuscript draft.Yong Pang and Wen Jia contributed to the data collection/validation and provided critical comments for the manuscript. Zehua Huang helped with data processing and provided some helpful comments for the manuscript. Feng Zhao conceived and designed the experiments, contributed to the data collection, allocated funding for the project,contributed to the first draft and revisions of the manuscript, and approved the final draft. All authors checked and approved the final manuscript.

Declaration of competing interest

The authors declare that they have no competing interests.

Acknowledgements

We would like to thank Zhengang Lv,Binyuan Xu,Yutao Zhao,Yigui Liao from Huazhong Agricultural University for their help during the field work. We also thank the staff members and managers from the Mengjiagang Forest Plantation for their local assistance.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://do i.org/10.1016/j.fecs.2022.100032.