A new strategy for improving the accuracy of forest aboveground biomass estimates in an alpine region based on multi-source remote sensing

ABSTRACT

Spatially explicit information on the distribution of dominant tree species groups and aboveground biomass (AGB) in forested areas is essential for developing targeted forest management and biodiversity conservation measures, as well as assessing forest carbon sequestration capacity. There is a shortage of continuously updated 30-m spatial resolution products for mapping dominant tree species groups. The vast majority of remote sensing-based AGB estimation approaches have relatively low accuracy for dominant tree species groups or forest types and are unsuitable for AGB modeling. Therefore, this study aims to develop an integrated framework that considers the phenological characteristics of different tree species to improve the mapping accuracies of forest dominant tree groups and corresponding AGB estimates. Thirty-meter resolution maps of dominant tree species groups were created using machine learning algorithms and phenological parameters. Features extracted from optical and radar images and phenological characteristics were used to construct AGB estimation models in a temporally consistent manner to improve the AGB estimation accuracy and perform dynamic AGB monitoring. The proposed method accurately characterized the dynamic distribution of the dominant tree species groups in the study area. The traditional AGB model that does not consider different forest types or species had an R² value of 0.52, whereas the proposed model that considers phenology and forest types had an R² value of 0.67. This result indicates that incorporating information on phenology and dominant species improves the accuracy of AGB estimations. The AGB in most regions was 30–55 t/ha, showing that the majority of the forests were young or middle-aged stands, and the areal percentage of AGB greater than 30 t/ha increased during the study period, suggesting an improvement in forest quality. Furthermore, the oak AGB was the highest, indicating that oak afforestation should be encouraged to enhance the carbon sequestration capacity of future forest ecosystems. The results provide new insights for researchers and managers to understand the trends of forest development and forest health, as well as technical information and a database for formulating more rational forest management strategies.

KEYWORDS:

• AGB
• forest species groups
• phenology
• remote sensing

1 Introduction

Spatially explicit information on the distribution of dominant tree species groups in forested areas is critical for numerous scientific applications, such as biodiversity protection, biomass estimation, and tree species management (Boschetti et al. 2007; Fassnacht et al. 2016; Shang and Chisholm 2014). Extensive field surveys conducted by forestry professionals provide highly accurate information. However, this approach is time-consuming, labor-intensive, and cost-ineffective, and areas with steep slopes can often not be accessed (Cao 2020). The use of modern remote sensing techniques coupled with well-tuned algorithms to derive dominant tree species distributions is a reliable alternative because of the wide availability of multi- and hyper-spectral, multi-temporal, and multi-scale resolution imagery. These approaches facilitate the extraction of forest information primarily based on spectral variables and textural features (Linderman et al. 2004; Zhu and Liu 2014; Araza et al. 2022; Fassnacht et al. 2016). Due to the high temporal revisit frequency of satellites such as Landsat, long-term observations have been widely used for large-scale dominant tree species classifications at 30 m spatial resolution. However, most studies have only distinguished broad classes of forest types, e.g. coniferous and broad-leaved forests, or have extracted single tree species (Linderman et al. 2004; Yang, Xu, and Zhai 2021). Detailed mapping of dominant tree species is generally performed using high spatial resolution or hyperspectral data (Cho, Malahlela, and Ramoelo 2015; Lin et al. 2015). Thus, there are insufficient mapping products for dominant tree species groups at medium spatial resolution.

The appropriate selection of variables or image features for dominant tree species group classifications is critical to improving the classification accuracy (Lu and Weng 2007). Phenological parameters are a promising feature for dominant tree species classification because the phenology of dominant tree species changes with the seasons. Some phenological events, such as leaves turning green in spring and yellow in autumn, and the timing of these events, are generally species-specific. Thus, tree species exhibit different spectral characteristics in multi-temporal remote sensing images, enabling the differentiation of dominant tree species in forestlands (Fassnacht et al. 2016; Yu et al. 2004). For example, Yu et al. (2004) used MODIS time-series data to classify forest types and found that phenological parameters improved the separability of forest types. Yang, Xu, and Zhai (2021) also confirmed that incorporating phenology into random forest (RF) modeling increased the mapping accuracy of dominant tree species. In addition to phenology, different spectral characteristics in different periods can be exploited to improve classification accuracy (Fassnacht et al. 2016). It has been suggested to avoid using late autumn images because the forest background signals make the separation of broadleaf species more difficult (Voss and Sugumaran 2008). Other studies have shown that the combination of images acquired in different seasons (one in spring and the other in autumn) can result in higher classification accuracy (Hill et al. 2010). In summary, we should consider phenological characteristics and the image acquisition time to classify dominant tree species.

Forest aboveground biomass (AGB) is critical to the global biogeochemical cycle and biodiversity conservation (Shen et al. 2016). Field measurement is the most accurate method for estimating forest AGB; however, this method lacks sufficient temporal and spatial coverage (Huiyi et al. 2020). As an alternative, Landsat imagery can be used to estimate forest AGB in large areas due to its adequate spatial and spectral resolution and wide coverage. In recent years, Landsat’s free and open policy has made long-term AGB assessments very popular. Many studies have shown that forest AGB and Landsat spectral data are highly correlated (Shen et al. 2016; Li, Li, and Li 2019; Huiyi et al. 2020). However, existing remote sensing-based AGB estimation models often have low accuracy due to an inadequate characterization of the forest structure. The distribution of tree groups is one of multiple horizontal stand structure characteristics and affects the accuracy of AGB modeling (Li et al. 2020; Li, Li, and Li 2019). Due to the lack of medium-resolution products suitable for characterizing dominant tree species groups in a changing landscape, most studies used a unified model for AGB estimation, regardless of the dominant tree species group or forest type, resulting in the overestimation (underestimation) of sparse (dense) forests (Singh et al. 2014). In addition, the selection of model variables is also highly critical for improving model accuracy. The most commonly used variables are spectral features extracted from optical remote sensing images, such as vegetation indices, and texture measures (Zhang et al. 2020; Bao et al. 2022; Pötzschner et al. 2022). However, these variables are limited because they lack temporal scalability and do not reflect the changing physiological characteristics of forests. New features are expected to improve AGB estimation. Research has demonstrated correlations between phenological characteristics and AGB (Ganjurjav et al. 2021; Keenan et al. 2014). For example, global warming prolongs the growing season and further increases forest carbon storage (Keenan et al. 2014). However, Ganjurjav et al. (2021) also found that the shortening of the growing season increased the biomass in an alpine meadow on the Tibetan Plateau. Therefore, capturing these seasonal changes via remote sensing can reveal physical changes in forests.

Developing reliable methods to extract phenological characteristics and dominant tree species groups is a promising option to improve the accuracy of AGB estimation models. The major contribution of this work to forestry remote sensing is the proposed progressive and time-consistent strategy for forest AGB estimation. Phenological characteristics are extracted to map the distribution of dominant tree species groups with high accuracy. Subsequently, features from optical and radar images are combined with the phenological features to construct a time-consistent AGB estimation model based on tree species groups instead of having to model AGB at multiple time points.

2. Materials and methods

2.1 Study area

The study area is located in the east of the Dali Bai Autonomous Prefecture (99°57 ”−101 °03‘ E, 25°14 ’−26 °13” N) of Yunnan Province, southwest China. It includes Dali City, Xiangyun County, Binchuan County, parts of Midu, Eryuan, and Heqing Counties, and the Yangbi and Weishan minority autonomous counties. The total area is 11,405 km² (Figure 1). The highest altitude in the study area is 4295.8 m, and the lowest is 730 m. The dry and wet seasons are distinct in this study area. Dry-season rainfall (from November to early April) accounts for 5–15% of the annual rainfall, whereas wet-season rainfall (from May to October) accounts for 85–95% of the annual rainfall. The average annual precipitation is approximately 1000 mm, and the annual average temperature is around 15°C (Ding and Peng 2018). Pines (Pinus yunnanensis, Pinus armandii Franch.), oak (Quercus pannosa), and walnut (Juglans sigillata Dode.) are the dominant tree species in the study area (Figure 2). The pine species are dominant in the evergreen coniferous forests, which are the most widely distributed forest type in this study area. These pine species have strong adaptability and can withstand cold, drought, and barren soil. Most oaks grow in evergreen broad-leaved forests and are adapted to the cold temperature at high altitudes. Walnut trees are a deciduous species with high economic value. They like sunlight and are suitable for rich and moist sandy soil. Although Dali prefecture has numerous vegetation species, the quality of some tree species is low, such as thin Pinus yunnanensis trunks and low density of walnut trees (Figure 2). Existing timber forests include natural and middle-aged and young plantation forests. Due to an imbalance between wood supply and demand, it is difficult to support the large-scale development of the forestry industry.

Figure 1. Study area; (a-c) show the geographical location of the study area and the phenological observation sites; (d) is a Landsat 8 OLI false-color composite (R: NIR band, G: red band, B: green band) acquired on 15 December 2019.

Figure 2. Field photos of common dominant tree species groups in the study area.

2.2 Datasets

We used Landsat and Sentinel-2 imagery, Phased Array L-band SAR (PALSAR) data, and forest field survey datasets. The Landsat surface reflectance (SR) data included blue, green, red, near-infrared (NIR), shortwave infrared 1 (SWIR1), and SWIR2 bands with less than 80% cloud cover ranging from 1999 to 2019. They were principally used to extract forest phenological parameters, including the start of the growing season (SOS), end of the growing season (EOS), and length of the growing season (LOS). These data were also used to map dominant tree species groups and AGB. Sentinel-2 top-of-atmosphere (TOA) reflectance images from 2017 to 2019 were used as supplementary data for Landsat imagery to increase the observation frequency and improve the extraction accuracy of the forest phenological parameters. The Sentinel-2 bands corresponded to those of Landsat. A total of 610 scenes were collected, including 514 scenes from Landsat Path/Row number: 131/042 (152 Landsat 5 TM images, 243 Landsat 7 ETM+ images, and 119 Landsat 8 OLI images) and 96 scenes from Sentinel-2 orbital number: T47RPJ (58 Sentinel-2A images and 38 Sentinel-2B images).

Five PALSAR images from 2007, 2010, 2015, 2016, and 2019 were used to extract features for the model input to characterize the vertical structure of the forest canopy and improve the AGB prediction accuracy.

Forest survey data were employed for the training and validation of the classifications of the dominant tree species groups and the forest AGB estimation. We collected data from the national forest inventory (NFI) plots in 2011 and forest management planning inventory (FMPI) in 2007 and 2016 (Lei et al. 2009; Xie et al. 2011). The NFI records various forest attributes (tree height, diameter at breast height, stand volume, canopy density, dominant species, etc.) using systematic sampling. The permanent plots with a size of 28.28 m × 28.28 m are revisited by local forest investigators every 5 years in China. A total of 78 permanent sample plots in the study area were measured in 2011. The FMPI divides the forest areas into multiple irregular sub-compartments that cover the entire research area. Local forest investigators record the average attributes of each sub-compartment nearly every ten years (Lei et al. 2009; Xie et al. 2011). A total of 2000 forest samples (0.1–2 ha) were randomly selected from the 2007 and 2016 data. It is worth noting that pure forest areas were selected as samples. CongDe et al. (2008) established a regression model between AGB and stock volume in Sichuan Province, southwest China, using extensive field measurements of Pinus yunnanensis, other pines, and oak. Thus, the species-specific biomass conversion factor was used to convert the stocking volume into forest AGB due to the geographic proximity (Appendix I) (CongDe et al. 2008) For dominant tree species without existing allometric biomass equations, the biomass equations of similar dominant tree species (same genus) were used to approximate their biomass. The field biomass data from 2007, 2011, and 2016 were combined to facilitate the establishment of the time-consistent estimation models.

2.3 Method

In our previous study (Zhang and Mingshi 2021), we integrated Landsat and Sentinel-2 to develop the modified continuous change detection and classification (MCCDC) model to extract phenological parameters. In this study, a similar method was again used to extract phenological parameters to further characterize the spatial distribution of forest tree species groups and AGB. Specifically, a series of multi-source data integration techniques were performed to match the Landsat and Sentinel 2 images after atmospheric correction and cloud mask (Zhang and Mingshi 2021). Subsequently, MCCDC model was adopted to synthesize cloud-free daily images and extract the MCCDC coefficients. Logistic regression models were used to extract the phenological parameters from the daily Landsat images (Zhang and Mingshi 2021). Based on MCCDC coefficients and phenological parameters, a map of dominant tree species groups was generated using a support vector machine (SVM). After the appropriate important variables were selected, a model with optical and PALSAR image features and phenological characteristics was constructed for different dominant tree species groups to estimate AGB in the entire study area. Figure 3 shows the flowchart of this study.

Figure 3. Flowchart of this study.

2.3.1 Image pre-processing

2.3.1.1 Pre-processing of Landsat and sentinel-2 data

Due to the large data volume, the atmospherically corrected Landsat SR data were ordered from the United States Geological Survey (USGS) Earth Resources Observation and Science Center (https://espa.cr.usgs.gov/ordering/new). The SR was multiplied by 10,000 to convert the data into an integer format for easy storage and handling. The USGS uses the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm for the atmospheric correction of TM and ETM+ images (Masek et al. 2006), whereas the Landsat 8 Surface Reflectance Code (LaSRC) algorithm is used for the atmospheric correction of the OLI images (Vermote et al. 2016). These two techniques are very mature and routinely implemented to generate highly accurate and consistent long-term (from 1982 to present) Landsat Analysis Ready Data (ARD) for surface change monitoring (Dwyer et al. 2018).

In addition, Sentinel-2 TOA reflectance data with less than 80% cloud cover from March to November 2017–2019 were obtained from the European Space Agency (ESA) as a supplement to Landsat data to increase the observation frequency. The bands were spectrally consistent with those of Landsat. The Sentinel-2 TOA reflectance data were atmospherically corrected to SR data using the Sen2cor radiative transfer atmospheric correction code (version 2.8, released 21 February 2019; http://step.esa.int/main/third-party-plugins-2/sen2cor/). The function of mask (Fmask) 4.0 software was used for batch cloud removal from all Landsat and Sentinel-2 data (Zhu and Woodcock 2014).

In addition to atmospheric correction and cloud masking, we also performed registration, resampling, BRDF correction, terrain correction, and bandpass adjustment to integrate the Landsat and Sentinel-2 data. Finally, Sentinel-2 was spatially and spectrally consistent with Landsat. The pre-processing parameters and results have been described by Zhang and Mingshi (2021).

2.3.1.2 PALSAR data pre-processing

The PALSAR images with a spatial resolution of 25 m acquired in 2007, 2010, 2015, 2016, and 2019 were used for the AGB estimation. The product contains horizontal-transmit, vertical receive (HV), and horizontal transmit, horizontal receive (HH) polarization images that have been geometrically and radiometrically corrected and terrain-normalized. The digital number (DN) was converted into the dB backscattering coefficients of the HH and HV polarization components as follows:

(1) ${\sigma }^0\left[dB\right]=10\times lo{g}_{10}\left(D{N}^2\right)+CF$(1)

where CF is the calibration factor. The CF for HH and HV were −83.2 and −80.2, respectively, before 2009 and −83.0 after 2009 (Shimada et al. 2009). An enhanced Lee filter with a window size of 5 × 5 was used to reduce the speckle noise (Jong-Sen et al. 2009). For most SAR data, the pixels in the 3 × 3 window may not be statistically independent due to oversampling to maintain the spatial resolution (Jong-Sen et al. 2009). Therefore, a 5 × 5 window size was used here. The related variables, including the HH and HV polarization values and the ratio and difference between HH and HV, were derived. All images were resampled to 30 m resolution using the nearest neighbor sampling method to be consistent with the Landsat images. All steps were implemented in the professional software The Environment for Visualizing Images (ENVI).

2.3.2 Extraction of forest phenological parameters

Zhang and Mingshi (2021) developed a new method to extract the SOS from multi-source remote sensing data (Landsat and Sentinel-2) using the MCCDC model. The MCCDC produces synthetic daily cloud-free Landsat images. The MCCDC model assumes that the SR value of pixels in areas with a stable land cover changes periodically over time. Thus, the MCCDC model first determines when a land cover change has occurred in a time series. It then simulates the periodic changes before and after the land cover changes and creates daily remote sensing images from the fitted curves. The synthesized SR images do not contain clouds, shadows, and scan lines.

Zhang and Mingshi (2021) also found that the land surface water index (LSWI) had a highest accuracy for extracting the SOS in this study area. Thus, a logistic regression model was used to extract the annual phenological parameters (SOS, EOS, LOS) using the LSWI in this study. When the rate of change of curvature of the logistic regression model reaches the first peak, the corresponding date (day of year, DOY) is the SOS, and when the rate of change of curvature reaches the last valley, the corresponding date is the EOS. The difference between the EOS and SOS is the LOS (Zhang et al., 2003). The MCCDC model was established, and the phenological parameters were extracted using the Interactive Data Language software.

2.3.3 Forest dominant tree group mapping based on the MCCDC coefficient and phenological parameters

2.3.3.1 Principles of forest dominant tree group mapping

We used cloud-free Landsat data (30 m resolution) to classify only the dominant tree species groups rather than specific tree species. After interpreting the 2016 FMPI data, we found that the dominant tree species in the study area were pine (53.52%), oak (30.20%), and walnut (8.5%). Thus, the forest was classified into four categories: pine, oak, walnut, and other dominant tree species. The spectral characteristics differed for different dominant tree species and were characterized by the MCCDC model coefficients (Figure 4). Specifically, the NIR reflectance of pine was substantially different from that of the oak and walnut species (Figure 4). The maximum SR of pine was generally lower than 3000, but the maximum SR of oak and walnut generally exceeded 3000. The spectral curves of walnut were different from those of the other two species in the red, SWIR1, and SWIR2 bands. Phenology can be used for dominant tree species classification because different tree species have different phenological characteristics. As shown in Figure 5, walnut has an earlier SOS than pine and oak, pine has a later EOS than oak and walnut, and walnut has a longer LOS than pine and oak. We divided the study area into forest and non-forest to facilitate the classification of the dominant tree species groups. We fed the regions of interest (ROIs) and MCCDC coefficients into an SVM to classify forest and non-forest areas (Burges 1998). Subsequently, the forest masks from the previous step were further subdivided into different dominant tree species groups. In this study, the ROIs of the dominant tree species group samples were randomly selected from the FMPI in 2007 and 2016. Due to the lack of FMPI data in the other years, persistent forest regions were extracted from the annual forest cover map from 1999 to 2019. We assumed no change in the dominant tree species groups in these persistent forest regions. The spatial intersections of the persistent forest areas and the samples from 2007 and 2016 were used to extract ROIs for the remaining years.

Figure 4. Spectral characteristics of different dominant tree species over time characterized by the MCCDC coefficients.

Figure 5. Comparison of the phenological parameters of different dominant tree species groups.

The MCCDC coefficients and phenological parameters were used as input into the SVM to obtain a map of the dominant tree species groups. Finally, the classification of the dominant tree species groups was performed in the kernlab program of RStudio.

2.3.3.2 Accuracy assessment of dominant tree species groups

Here, a ten-fold cross-validation strategy was used to evaluate the classification accuracy of the model (Fushiki 2011). All ROIs were divided into 10 parts; nine parts were used as training data, and one part was used as validation data. This process was repeated 10 times, and majority voting was implemented to generate the final classification result and calculate the accuracy statistics of the algorithm. The accuracy was evaluated using the user’s accuracy, producer’s accuracy, and overall accuracy. The user’s accuracy refers to the ratio of correctly classified samples to the classified samples of a certain category. The producer’s accuracy represents the ratio of correctly classified samples to the number of reference samples of a certain category. The overall accuracy refers to the ratio of correctly classified samples to the samples of all categories. Additionally, the binomial proportion confidence interval, which is based on the Bernoulli trial, was used to express the uncertainty of the classification accuracy (Eq. 2) (Orawo 2021).

(2) $d=P\pm Z\times \sqrt{P\left(1-P\right)/n}$(2)

Where $\mathrm{d}$ is the confidence interval, P refers to the estimated accuracies, including the user’s accuracy, producer’s accuracy, and overall accuracy, n is the number of validation samples. Z equals 1.96 when the confidence level is 95%.

2.3.4 Establishment of forest AGB estimation model

2.3.4.1 AGB prediction variables

The model variables included optical remote sensing and topographic variables, including SR, normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), soil-adjusted vegetation index (SAVI), tasseled cap transformation, principal component analysis (PCA), texture information, elevation, slope, and aspect (Byrne, Crapper, and Mayo 1980; Crist, Eric, and Richard 1984; Haralick, Shanmugam, and Dinstein 1973; Huete et al. 2002; Huete 1988; Schell 1973). The tasseled cap transformation creates a new set of features using a linear combination of six Landsat SR bands and fixed coefficients. The first three features are brightness, wetness, and greenness, which have physical meanings in vegetation analysis (Crist, Eric, and Richard 1984). PCA is an orthogonal linear transformation to remove redundant information and recombine the variance contributions into newly generated components (Byrne, Crapper, and Mayo 1980). We also used the gray-level co-occurrence matrix and six SR bands to calculate the mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, and correlation (Haralick, Shanmugam, and Dinstein 1973). These operations were implemented in ENVI. Elevation data with a spatial resolution of 90 m were obtained from Shuttle Radar Topography Mission data from the geospatial data cloud website (http://www.gscloud.cn/). The elevation data were resampled to 30 m spatial resolution using the nearest neighbor resampling method. The slope and aspect tools of ArcGIS were used to obtain the slope and aspect data layers. We incorporated the MCCDC coefficients and phenological parameters (SOS, EOS, LOS) into the AGB model to improve the AGB estimation accuracy. The PALSAR data and its derivatives (HH, HV, and their difference and ratio) were also included in the AGB model (Table 1).

Table 1. Predictor variables for AGB modeling.

Type	Variables
SR data	Blue, red, green, NIR, SWIR1, SWIR2
Vegetation indices	NDVI
	EVI
	SAVI
Tasseled cap transformations	Brightness, wetness, greenness
Texture measures	Mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, correlation
PCA	PCA variables
Topographic variables	Elevation, slope, aspect
Seasonal variation	MCCDC coefficients
Phenological parameters	SOS, EOS, LOS
PALSAR data	HH, HV, and the difference and ratio of HH and HV

2.3.4.2 Selection of important variables and AGB modeling

Using redundant variables in an AGB model reduces its prediction ability (Shen et al. 2016). Therefore, the variable importance function of the RF package was used to rank the importance of the variables using two indicators: the percent increase in the mean square error (PercentIncMSE) and the increase in node purity (IncNodePurity) (Breiman 2001). The higher the PercentIncMSE and IncNodePurity are, the more important the variables are.

The important predictor variables and the matching measured AGB data were used to establish the time-consistent AGB prediction models at three-year interval in the RF package of R software. Three parameters were set in this model: 1) the ntree is the number of decision trees; it was set to 500; 2) the mtry is the number of variables at each node; it was set to one-third of the input variables; 3) the node size is the minimum number of nodes in the decision tree; we used the default value of 1.

2.3.4.3 Evaluation of forest AGB model

Similar to the evaluation of the dominant tree species groups classification, we used ten-fold cross-validation to evaluate the prediction accuracy of the AGB prediction model (Fushiki 2011). The results from ten outputs were averaged. The R², root mean square error (RMSE), mean absolute error (MAE), and RMSE percentage (RMSE%) between the actual and predicted values were selected to evaluate the model performance (Eq. 3–6). In general, the larger the R², and the smaller the RMSE, MAE, and RMSE%, the better the predictive ability of the model.

(3) $R^2=1-\frac{{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y_i}\right)}^2}$(3)

(4) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\sum }_{i=1}^n\frac{{\left(y_i-\widehat{y_i}\right)}^2}{n}}$(4)

(5) $\mathrm{M}\mathrm{A}\mathrm{E}=\frac{1}{n}{\sum }_{i=1}^n\left|\widehat{y_i}-y_i\right|$(5)

(6) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}\mathrm{\%}=\frac{RMSE}{\bar{y}}\times 100$(6)

where $y_i$ is the measured AGB, $\widehat{y_i}$ is the estimated AGB, $\bar{y}$is the mean of the observed AGB values, and n is the number of samples.

3 Results

3.1 Classification accuracy and distribution of forest dominant tree species groups

The accuracy of the forest phenology estimation is shown in Appendix II. Table 2 and Table 3 indicate that the classification accuracy of the model with the MCCDC coefficients and phenological parameters was higher than that of the spectral classification algorithm that used SR, texture information, EVI, and NDVI. A t-test was used to determine the statistically significant difference in the model prediction capability of the two methods. The differences between the left and right endpoints of the confidence intervals did not exceed 3% for all accuracies (Table 2 and Table 3). Thus, the MCCDC coefficients and phenological parameters were used as the optimal input variables to create the maps of the forest dominant tree species groups in the study area.

Table 2. Accuracy and confidence interval of the dominant tree species groups obtained from the model with only spectral variables in 2007 and 2016 (confidence level = 95%).

	2007		2016
	User’s accuracy	Producer’s accuracy	User’s accuracy	Producer’s accuracy
pine	75%±1.2%	85%±1.0%	74%±1.2%	80%±1.1%
oak	79%±1.4%	80%±1.4%	74%±1.5%	74%±1.5%
Walnut	71%±2.3%	43%±2.5%	64%±2.3%	45%±2.0%
other dominant tree species	84%±1.7%	28%±2.1%	78%+±2.5%	34%±2.9%
non-forest	89%±1.0%	91%±0.9%	88%±1.0%	91%±0.9%
overall accuracy	80%±0.7%		78%±0.7%

Table 3. Accuracy and confidence interval of the dominant tree species groups obtained from the model with the MCCDC coefficients and phenological parameters in 2007 and 2016 (confidence level = 95%).

	2007		2016
	User’s accuracy	Producer’s accuracy	User’s accuracy	Producer’s accuracy
pine	87%±1.4%	92%±0.8%	87%±1.0%	88%±1.1%
oak	88%±1.1%	91%±1.0%	82%±1.3%	93%±0.9%
walnut	91%±1.4%	86%±1.7%	91%±1.4%	87%±1.6%
other dominant tree species	87%±1.6%	74%±2.0%	88%±2.0%	71%±2.8%
non-forest	94%±0.7%	91%±0.9%	94%±0.7%	90%±0.9%
overall accuracy	89%±0.5%*		88%±0.5%*

*indicates significant differences in the accuracy of the two methods at the 95% confidence level.

Figure 6 shows the spatial distribution of the dominant tree species groups from 2000 to 2019. Pines were more common on the right side of the western mountains and the edge of the other mountains. Oaks occurred predominantly on the left side of the western mountain and the upper right and middle upper part of the study area. Walnuts were mainly located at the left edge of the study area. In the upper right corner of the study area, oaks were interspersed with walnuts, and other forest species were scattered at the edges and inside other forest stands.

Figure 6. Spatiotemporal distribution of the dominant tree species groups obtained from the model with the MCCDC coefficients and phenological parameters.

Figure 7 shows the percentage of different dominant tree species groups in different years. The area of pine forest decreased over time, except for a temp orary increase from 2007 to 2010. However, the trend of the oak forest was the opposite of that of the pine forest. The percentage of walnut forest remained unchanged at about 12%, and the percentage of other dominant tree species groups showed a downward trend.

Figure 7. Percentage of different dominant tree species groups from 2000 to 2019.

3.2 Forest AGB validation accuracy

3.2.1 Effects of different variable combinations on AGB estimation accuracy

A total of 16 variables were selected for this study according to the order of the two indicators (Red_3, NIR, NIR_1, NIR_3, SWIR1, SWIR1_3, SWIR2_5, PCA1, PCA3, HV, HHHV, SOS, EOS, LOS, DEM, Wetness) (Figure 8). Note that NIR_1 refers to the first coefficient of the MCCDC for the NIR band, and so on.

Figure 8. Variable importance ranking of the predictor variables.

Different variable combinations were evaluated to demonstrate the influence of the MCCDC coefficients and phenological parameters on the AGB model. The RF model with only the spectral and PALSAR variables obtained an R² of 0.52 (Table 4). In contrast, the model that included the MCCDC coefficients and phenological parameters had an R² 0.56. The RMSE was 18.62 t/ha for the first model and 14.05 t/ha for the proposed model. The MAE was slightly higher for the proposed model, and the RMSE% was 52.49% for the first model and 32.91% for the proposed model (Table 4). Therefore, incorporating the MCCDC coefficients and phenological parameters resulted in a higher accuracy of the AGB model.

Table 4. Evaluation indices of models with different variable combinations.

	R^2	RMSE (t/ha)	MAE (t/ha)	RMSE%
Spectral variables+PALSAR	0.52	18.62	13.66	52.49
Spectral variables + phenology variables+PALSAR	0.56	14.05	13.12	32.91

3.2.2 Influence of dominant tree species on AGB estimation accuracy

The results revealed the validation accuracy of the AGB model when separating tree species groups for AGB modeling. Walnut had the highest validation R² (0.78), with relatively small RMSE, MAE, and RMSE%. The R² values of the pine, oak, and other dominant tree species were 0.61, 0.62, and 0.65, respectively. Although the R² of the other dominant tree species was slightly higher than those of pine and oak, the RMSE, MAE, and RMSE% of this group were the highest (24.45 t/ha, 18.62 t/ha, and 56.67%, respectively) (Table 5).

Table 5. AGB accuracy of the proposed AGB model for different dominant tree species.

	R^2	RMSE（t/ha）	MAE（t/ha）	RMSE%
pine	0.61	9.75	8.35	32.47
oak	0.62	12.92	16.43	25.99
walnut	0.78	4.80	3.69	35.13
other dominant tree species	0.65	24.45	18.62	56.67

Figure 9 presents the scatter plots of the observed and predicted AGB for all samples before and after considering the dominant tree species groups. There were significant differences in the predicted AGB between the two methods (Figure 9). The R² and RMSE were 0.56 and 14.05 t/ha in the former case and 0.67 and 12.13 t/ha in the latter case. Some low (high) AGB values were overestimated (underestimated) prior to considering the tree species, particularly when the AGB <15 t/ha and >60 t/ha. However, the differences were less pronounced when the dominant tree species were considered because the y=x line of the AGB model was closer to the trend line. This result demonstrates that each dominant tree or forest type is required for AGB model estimation.

Figure 9. Scatter plot of predicted and observed AGB; (a) without considering dominant tree species and (b) considering dominant tree species. * indicates significant differences in the predicted AGB between (a) and (b).

3.3 Temporal and spatial dynamics of AGB

We used the proposed models for different dominant tree species to map the spatial distribution of AGB from 2000 to 2019. The AGB in the lower left region of the study was less than 30 t/ha, and walnut was the dominant tree species in this area (Figures 6 and 10) In the southeast corner of the study area, the AGB was generally lower than 30 t/ha, and pine was dominant. The AGB was relatively high (greater than 55 t/ha) on the left side, the northeast corner, and the upper middle area of the study area, and oak was the main species. In the northwest corner and the lower central region, the AGB increased from 30–55 t/ha to more than 55 t/ha from 2000 to 2019. The AGB in the southeast corner did not change substantially over time and was less than 30 t/ha (Figure 10).

Figure 10. Spatial distribution of forest AGB from 2000 to 2019.

Figure 11 presents the AGB statistics for different years. From 2000 to 2019, the area percentage of AGB less than 30 t/ha exhibited a downward trend. From 2000 to 2013, the area percentage of AGB of 30–55 t/ha increased, reaching 45.83%. The area percentage of AGB of 55–80 t/ha increased from 12.54% in 2000 to 18.40% in 2019, and the area percentage of AGB of 80–105 t/ha increased from 4.19% in 2007 to 9.76% in 2016. The area percentage of AGB greater than 105 t/ha remained less than 0.8%. The area percentage of AGB of 30–105 t/ha increased, whereas that with an AGB of less than 30 t/ha decreased. These results indicate an increase in forest quality in the study area from 2000 to 2019.

Figure 11. AGB statistics from 2000 to 2019.

Figure 12 shows the total AGB of different dominant tree species in different years. The total AGB of pines remained unchanged from 2000 to 2007 and decreased after 2010 to 9.57E + 06 ton in 2019. The total AGB of oak increased to 1.20E + 07 ton, except for a decline in 2010. The total AGB of walnut was lower than that of pine and oak. The total AGB of other forest species decreased from 2.93E + 06 ton to 2.39E + 06 ton from 2000 to 2007 and then remained at 2.90E + 06 ton.

Figure 12. Total AGB of different dominant tree species from 2000 to 2019.

4 Discussion

4.1 Effects of phenological characteristics on dominant tree species group classification

The phenological characteristics differed for different species. The SOS was generally later for pine and oak than for walnut, pine had a later EOS than oak and walnut, and walnut had a longer LOS than the other dominant tree species (Figure 5). Besides, different dominant tree species exhibited different spectral response curves over time. Therefore, incorporating the MCCDC coefficients and phenological parameters into the model improved the classification accuracy of the dominant tree species groups.

As shown in Figures 6 and 7, pines were the most common tree species in the study area, followed by oaks, although the area of the oak forest increased over time. Pinus yunnanensis is a native dominant tree species in southwest China, with good adaptability and drought resistance (Cha 2018); therefore, pines were the most widely distributed tree species. However, a single-species coniferous forest does not have high biodiversity and is easily threatened by forest insects and disease; thus, it is necessary to increase the proportion of broad-leaved trees to create more stable mixed forest ecosystems. The afforestation project has increased the proportion of oaks in recent years (https://www.docin.com/p-1263061185.html).

4.2 Advantages of estimating forest AGB based on dominant tree species and seasonal characteristics

Our results showed that the validation R² of the AGB model was higher (0.56 compared to 0.52 for the other model) after the phenological parameters were added. Some studies have shown that global warming advances the SOS and delays the EOS, lengthening the LOS and increasing AGB (Chu et al. 2016; Jin et al. 2019). However, a delay in the SOS in our study area also increased AGB. The reason may be that the growing season was shortened, which decreased the drought duration, leading to an increase in AGB (Figures 11 and 12)

Different dominant tree species have different AGB values. Our results indicated that oak species had higher AGB, while walnut species had lower AGB, even though they are both broad-leaved tree species. In addition, pine had lower AGB than oak (Figures 6 and 10). Thus, separate AGB models for different forest types are better able to capture the differences in canopy characteristics between different coniferous and broadleaved forests. The validation of the AGB estimation accuracy showed that when one AGB model was used for all tree species, low AGB values (less than 15 t/ha) were overestimated, and high AGB values (more than 60 t/ha) were underestimated. The AGB model that considered different dominant tree species had a higher R² and a lower RMSE, and the trend line approached y=x (Figure 9). Thus, it is necessary to construct AGB models for different tree species.

4.3 Comparison and drivers of AGB dynamics

(Han 2014) estimated forest AGB using one-year Landsat TM images and found that most AGB estimates were in the range of 35–55 t/ha. In contrast, Huang et al. (2013) obtained AGB estimates of 0–25 t/ha in the same area in 2005. The AGB values in our study area were 30–55 t/ha, suggesting low AGB. In addition, Cha (2018) estimated that the average forest carbon stock in the Erhai basin was 39.63 t/ha and found that the order of carbon density was oak > pine > other broad-leaved trees, which is consistent with our results.

In general, the AGB in this study area was low (Figures 10 and 11), which can be attributed to the wide distribution of young and middle-aged forests, poor site conditions, and ineffective forest management practices. For example, the main tree species is Pinus yunnanensis, which is managed for timber production with a short rotation cycle. Thus, frequent logging and regeneration events result in a positively skewed distribution of forest age (Cha 2018). Additionally, the walnut AGB was also low. The trunks are small, and the tree crowns are pruned during cultivation. In addition, the row spacing of the walnut trees is wide, resulting in a low number of walnuts per unit area. These factors are responsible for the low walnut AGB. However, the area percentage of AGB greater than 55 t/ha for oak increased from 2000 to 2019 (Figure 11) because of afforestation (Figure 7), leading to an increase in overall forest AGB. Besides, some young oak forests shifted to middle-aged and mature forests, increasing AGB. This information can be used to learn from past operational practices and improve forest management plans to achieve sustainable development in the near future. For example, afforestation is critical for increasing forest carbon sinks. Since oak has the highest AGB per unit area, it is suggested to prioritize oak afforestation in the study area to change the forest structure, which is currently dominated by coniferous forests, and achieve a transformation from single-species to multi-species forests. Mature forests generally have a larger carbon sequestration capacity. Thus, it is necessary to promote forest succession, prevent reverse succession, and focus on cultivating mature forests. Lastly, the greater the forest volume per unit area, the greater its carbon sequestration capacity. The most effective strategy to increase forest stock per unit area is to breed more long-lived, large-diameter trees (Yan et al. 2016).

4.4. Research limitations and future work

First, the Sentinel-2 data were resampled to 30 m resolution to match the Landsat resolution. In this process, some non-forest pixels might be mixed with forest pixels at the edge of forest and non-forest areas, decreasing the purity of the Sentinel-2 pixels. Therefore, some mixed pixels may have resulted in spectral information that did not accurately reflect the seasonal changes of the tree species. Second, cloud-free observations are required to establish the MCCDC model. In regions with cloudy conditions in the growing season, there may be an insufficient number of cloud-free observation points, potentially reducing the fitting accuracy of the MCCDC model. Thus, these MCCDC coefficients will affect the accuracy of tree species classification and AGB estimation. For example, we obtained a relatively low producer’s accuracy of about 74% for the other dominant tree species (Table 3). There was a considerable amount of misclassification of this class. The current availability of Landsat 9 and the upcoming launch of Sentinel-2C will provide more cloud-free data, improving the applicability and accuracy of the proposed method. Third, some high AGB values (>60 t/ha) were underestimated, and low AGB values (<15 t/ha) were overestimated (Figure 9), which was caused by the inadequate characterization of the forest structure, particularly the vertical structure. Using forest canopy height data obtained from lidar data to model the vertical stand structure may minimize this problem.

5 Conclusion

Incorporating phenological parameters and dominant tree species groups is considered an effective method to improve the estimation accuracy of AGB models. Therefore, this study created maps of dominant tree species groups and established a model that included features extracted from optical and radar images and phenological characteristics to improve the accuracy of forest AGB estimation. The results indicated that incorporating phenological parameters and dominant tree species data significantly improved the accuracy of estimating forest AGB. The area percentage with AGB greater than 55 t/ha increased from 2000 to 2019, and the oak AGB was the highest in the study area. The findings indicate that the forest quality level is improving. Oak afforestation should be encouraged. Another contribution of this work is that the proposed time-consistent AGB modeling strategy can be used to interpolate or extrapolate AGB estimates.

Acknowledgments

Special thanks are extended to the USGS EROS Data Center and ESA for providing the Landsat and Sentinel-2 images used in this work.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper

Data availability statement

The data that support the findings of this study (i.e. the land cover and probability maps) are available from the corresponding author and USGS (https://glovis.usgs.gov/).

Additional information

Funding

This work was jointly funded by the National Natural Science Foundation of China under grant number 31971577, the Foundation of Anhui Province Key Laboratory of Physical Geographic Environment under grant number 2022PGE011, the General Items of Support Plan for Outstanding Young Talents in Colleges and Universities under grant number gxyq2021217, the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD), and the Open Fund of State Laboratory of Information Engineering in Surveying, Mapping, and Remote Sensing, Wuhan University, under grant number 21R04.

Appendix Appendix I

Relationship between biomass and tree volume (Huang et al. 2008)

Table

Dominant tree species	biomass conversion equation
Pinus yunnanensis	w = 0.8569 v^0.8564
Other pines and conifer forests	w = 0.5272 v^1.0793
Quercus and Hardwood	w = 0.6573 v^1.0502
Eucalyptus	w = 2.5878 v^0.6261
Populus	w = 0.7629 v^0.9505
Softwood	w = 2.1529 v^0.6085
Abies and Picea	w = 2.137 v^0.7532
Cunninghamia lanceolata	w = 0.1807 v^1.2771
Cupressus funebris	w = 1.2237 v^0.8592
Cryptomeria fortunei	w = 0.8707 v^1.0239
Phoebe and Cinnamomum	w = 1.0726 v^1.0128

Appendix II

Scatter plots of observed and predicted SOS and EOS extracted from the EVI, NDVI, and LSWI