Mapping Spatiotemporal Changes in Forest Type and Aboveground Biomass from Landsat Long-Term Time-Series Analysis—A Case Study from Yaoluoping National Nature Reserve, Anhui Province of Eastern China

Abstract

A natural reserve’s forest is an important base for promoting natural education, scientific research, biodiversity conservation and carbon accounting. Dynamic monitoring of the forest type and forest aboveground biomass (AGB) in a nature reserve is an important foundation for assessing the forest succession stage and trend. Based on the Landsat images covering the National Nature Reserve of Yaoluoping in Anhui province spanning from 1987 to 2020, a total of 42 Landsat scenes, the forest cover product set was first developed by using the well-established vegetation change tracker (VCT) model. On this basis, a new vegetation index, NDVI_DR, which considers the phenological characteristics of different forest types, was proposed to distinguish coniferous forest from broad-leaved forest. Next, multiple modeling factors, including remote sensing spectral signatures, vegetation indices, textural measures derived from gray level co-occurrence matrix and wavelet analysis and topographic attributes, were compiled to model the AGB in 2011 by forest type separately by using the stochastic gradient boosting (SGB) algorithm. Then, using the 2011 Landsat image as the base, all the Landsat images in the other years involved in the modelling were relatively normalized by using the weighted invariant pixels (WIP) method, followed by an extrapolation of the 2011 AGB model to other years to create a time-series of AGB. The results showed that the overall accuracy of the VCT-based forest classification products was over 90%. The annual forest type classifications derived from NDVI_DR thresholding gained an overall accuracy above 92%, with a kappa coefficient above 0.8. The 2011 forest-type-dependent SGB-based AGB estimation model achieved an independent validation R² at 0.63 and an RMSE at 11.18 t/ha for broad-leaved forest, and 0.61 and 14.26 t/ha for coniferous forest. The mapped time-series of AGB showed a gradual increasing trend over the past three decades. The driving factors responsible for the observed forest cover and AGB changes were analyzed to provide references for reasonable protection and development. The proposed methodology is a reliable tool for evaluating the management status, which can be extended to other similar regions.

1. Introduction

Forest provides 80% of the global aboveground vegetation biomass [1]. Obtaining reliable long-term forest change information over wide regions in an efficient, low-cost and timely manner becomes one of the major missions of the modern remote sensing framework to satisfy the technical demands of sustainable forest management. Particularly, persistent and accurate estimation of forest AGB or carbon storage in Chinese ecosystems has been a technical and data basis for achieving the “carbon peak before 2030 and carbon neutrality before 2060” aim [2,3,4].

Due to poor timeliness and limited accessibility, traditional forest field inventory has difficulty meeting the requirements of temporal and spatial dynamic monitoring of forest resources in a wide region and a long time span [5]. With the rapid development of modern remote sensing technology and computerized image analysis algorithms, multi-temporal and multi-resolution remote sensing images provide an important data source for land cover change studies at the landscape, regional and global scales, of which forest change information extraction based on long-term time-series remote sensing observations has attracted more and more attention of scholars [6,7,8,9,10]. The early study of forest change by remote sensing is generally based on the pairwise comparison of images between two periods [11,12,13]. However, this method may miss those fast-change events, for example, forest loss caused by fire or harvesting and forest gain due to immediate regeneration following harvesting in southern China, due to the relatively long time interval of the two periods (e.g., 10 years). Therefore, to adequately characterize forest changes, we need to use dense time-series images, for example, yearly observations. The Landsat imagery has a medium spatial resolution (30 m) and long-term archive of data (around 50 years), which is suitable for creating a long-term dense time-series stack to sufficiently monitor forest change at the landscape scale. Thanks to its free availability and long historical archive, Landsat-based long-term forest change analysis has been widely implemented by using several reputational automated analytical algorithms, including Landtrendr [8], VCT [9] and Breaks For Additive Season and Trend Monitor (BFAST) [10]. Among them, VCT has attracted much attention in the long-term forest change analysis because of its advantages of automation, high efficiency, high accuracy, easy implementation and full utilization of time information.

Mapping the forest type from remote sensing is one of the important contents of forest resources investigation and monitoring, which can help assess the successional stages and trends of a particular ecosystem. Additionally, the existing studies have shown that the accuracy of separately modeling the AGB by different forest types is generally higher than that of modeling the AGB without differentiating forest types or combining all the forest types together [14,15,16]. However, in practice, spatially explicit data of forest type are not always available or not suited to support the AGB modeling. For example, although the existing land cover products, e.g., the American NLCD and China’s land use datasets, generally contain forest types data (coniferous forest, broad-leaved forest and mixed forest classes), these products generally have a relatively long time cycle, every five years, to update; thus, they cannot adequately record forest type changes induced by frequent forest harvesting and regeneration events occurring within the time interval, let alone accurately capture the phenological differences of forest types [17,18]. Furthermore, these products generally have limitations on local-scale uses due to their inaccuracy at this scale [19]. At present, methods such as decision tree, support vector machine, random forest and neural network have been widely used in the field of forest type classification research to satisfy specific or personalized needs at the local or regional scale [20,21,22,23]. However, these methods generally require a large number of training samples to train the classification models, and the model’s parameters must be elaborately tuned to gain a high accuracy of classifications, which result in these methods being less efficient and not easily implemented. Thus, developing accurate and efficient methods to map forest types to facilitate the accurate modelling of dense time-series AGB, e.g., annual AGB, is of high priority.

Forest AGB is an underlying indicator for evaluating forest carbon sequestration capacity and biodiversity carrying capacity. The traditional field method for AGB estimation is subject to temporal and spatial limitations [24]. Currently, remote sensing combined with statistical modeling technology acts as a promising alternative to provide macro, near-real-time and multi-scale forest AGB estimation products. The earlier empirical modeling methods include simple linear regression or multiple linear regression models. These methods require a strict normal distribution assumption regarding the used dataset, but remote-sensing-derived variables and other GIS layers often do not conform to this distribution, so their use is greatly limited [25,26]. With the rapid development of computer technology, artificial intelligence has attracted more and more attention of scholars [27]. Machine learning, as one of the important components of artificial intelligence, with the advantages of not requiring data distribution assumptions, nonlinear complex mapping and being insensitive to sample outliers [28], with high model estimation accuracy through self-learning [29,30], has been predominantly used in classification and modelling applications. For example, Lawrence et al. [31] compared the classification results of the SGB algorithm and decision classification trees based on three kinds of datasets, indicating that the SGB algorithm can improve the classification accuracy. Guneralp et al. [32] used Landsat 7 ETM+ and SPOT 5 data to compare the accuracy of SGB-modeled forest AGB, multivariate adaptive regression spline algorithm and Cubist algorithm and indicated that the SGB algorithm was more accurate after adding terrain data, such as elevation, slope and aspect. Dube et al. [33] constructed texture and spectral features based on Landsat images and found that the SGB algorithm had higher accuracy in AGB estimation than random forest. Therefore, coupling the SGB algorithm and Landsat data can generate more accurate AGB inversion results than other machine learning algorithms and classic statistical methods. However, it is noted that these modelling methods are just used in very limited time points because adequate field forest measurement sample plots collected in multiple years, such as the model training set, are not available; thus, dynamically mapping AGB in a dense-time-series manner remains extremely difficult or even impossible. Therefore, this challenge necessitates the development of an effective and reliable framework that integrates long-term Landsat observations, the SGB algorithm and in situ biomass measurements in one year to dynamically model AGB.

Yaoluoping National Nature Reserve is a typical representative base for biodiversity conservation in the Dabie Mountains, Anhui province of eastern China. It is known as the “gene bank of natural species”. Understanding its changes in forest resources can help formulate more scientific and reasonable development and protection policies or actions. Therefore, the main objectives of this paper were to develop an efficient and reliable framework to generate multi-temporal AGB, and to investigate the driving factors responsible for forest temporal and spatial changes to recommend targeted policy suggestions for better management of the reserve. Specifically, the major contributions of the current work lie in: (1) developing an efficient and accurate image index, NDVI-DR, to map forest types in multiple years, and (2) devising a new framework that considers forest types, the SGB algorithm, Landsat time-series observations and one-year field biomass measurements to extrapolate the established 2011 biomass estimation model to other years in order to realize the dynamic generation of forest biomass products.

2. Materials and Methods

2.1. Study Area

The Yaoluoping National Nature Reserve, covering a total area of 123 km², is located in the northwest of Yuexi County, Anqing City, Anhui Province, with latitudes of 30°57′20″ to 31°06′10″N and longitudes of 116°02′20″ to 116°11′53″E (Figure 1). It has an average elevation of 800 m and belongs to the North subtropical monsoon zone. Because its location is between the Yangtze River and the Huaihe River, the reserve is affected by cyclones, with an abundant precipitation. According to the records of China Meteorological Administration, the annual mean temperature of the reserve is about 12.7℃ and the annual rainfall is 1700 mm. Yaoluoping was approved by the State Council to establish a national nature reserve on 5 April 1994. In 1999, Yaoluoping National Nature Reserve was included into the “Chinese people and Biosphere” list, with a dominant orientation of “Forest ecology”. More than 40 national key rare and endangered animal and plant species in the reserve had become important protection objects due to the inclusion of membership. The vegetation in the reserve is mainly deciduous broad-leaved forest and evergreen coniferous forest, and the dominant species are Cunninghamia lanceolata (lamb.) Hook., Pinus massoniana Lamb., Pinus dabeshanensis Cheng et Law, Pinus taiwanensis Hayata, Quercus stewardii Rehd. Quercus variabilis Bl., Alnus trabeculosa Hand.-Mazz. [34].

2.2. Data and Preprocessing

The remote sensing data used in this study included 42 scenes of Landsat TM/OLI imagery from 1987 to 2020, with a WRS-2 path/row number 122/038. The detailed description of the images was summarized in

Table 1. Description of the Landsat TM/OLI images used in the current work.

Acquisition Date	Satellite/Sensor	Cloud %	Acquisition Date	Satellite/Sensor	Cloud %
7 February 1987	Landsat 5 TM	0%	17 September 2004	Landsat 5 TM	0%
19 September 1987	Landsat 5 TM	0%	18 July 2005	Landsat 5 TM	35%
21 September 1988	Landsat 5 TM	9%	19 June 2006	Landsat 5 TM	1%
23 August 1989	Landsat 5 TM	33%	29 January 2007	Landsat 5 TM	0%
23 June 1990	Landsat 5 TM	40%	25 August 2007	Landsat 5 TM	23%
29 August 1991	Landsat 5 TM	1%	12 September 2008	Landsat 5 TM	49%
5 February 1992	Landsat 5 TM	9%	17 October 2009	Landsat 5 TM	0%
30 July 1992	Landsat 5 TM	1%	4 October 2010	Landsat 5 TM	5%
3 September 1993	Landsat 5 TM	8%	24 January 2011	Landsat 5 TM	2%
24 October 1994	Landsat 5 TM	0%	5 September 2011	Landsat 5 TM	69%
5 June 1995	Landsat 5 TM	1%	12 October 2013	Landsat 8 OLI	0.04%
25 July 1996	Landsat 5 TM	1%	25 January 2013	Landsat 8 OLI	2.1%
18 February 1997	Landsat 5 TM	0%	15 October 2014	Landsat 8 OLI	1.91%
29 August 1997	Landsat 5 TM	2%	2 October 2015	Landsat 8 OLI	0.01%
28 May 1998	Landsat 5 TM	23%	4 October 2016	Landsat 8 OLI	20.39%
2 July 1999	Landsat 5 TM	2%	25 February 2017	Landsat 8 OLI	11.84%
18 June 2000	Landsat 5 TM	0%	7 October 2017	Landsat 8 OLI	17.18%
23 July 2001	Landsat 5 TM	19%	8 September 2018	Landsat 8 OLI	0.77%
30 December 2001	Landsat 5 TM	2%	27 September 2019	Landsat 8 OLI	15.83%
28 September 2002	Landsat 5 TM	5%	18 February 2020	Landsat 8 OLI	2.12%
29 July 2003	Landsat 5 TM	16%	28 August 2020	Landsat 8 OLI	1.14%

. All images were downloaded from the USGS official portal (https://glovis.usgs.gov/, accessed on 20 January 2022). In order to explore the feasibility of using seasonal differences to distinguish coniferous forest from broad-leaved forest, images acquired in winter season were also downloaded for partial years of 1987, 1992, 1997, 2002, 2007, 2011, 2013, 2017 and 2020. Before placing the data order to USGS EROS, we directly requested the Landsat Level-2 products to lower the preprocessing complexity or workload. The Level-2 products are time-series observational data of sufficient length, consistency and continuity to record effects of climate change, and they are research-quality, applications-ready and generated for viable surface reflectance science data by USGS EROS data center. Specifically, Landsat 8 Operational Land Imager (OLI) surface reflectance products are generated using the Land Surface Reflectance Code (LaSRC) algorithm [35]. Landsat 5 TM surface reflectance products are generated using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [36].

Other auxiliary data used in this study included the 2012 vegetation distribution map of Yaoluoping National Nature Reserve, the 2011 forest resource inventory and planning data, the administrative boundary and functional zoning vector files, DEM data with a spatial resolution of 12.5 m and the Statistical Yearbooks of Anqing City in 2011 through 2019. Additionally, Google Earth maps were also collected to support the validation of classifications.

2.3. VCT-Based Forest Distribution Extraction in Yaoluoping Nature Reserve

VCT algorithm was developed by Huang et al. at the University of Maryland in 2010 [9], and it has been widely used and tested around the world in recent decades to characterize forest change patterns, with an average overall accuracy at about 85% [37,38,39]. Here, we directly applied this algorithm to create a time-series of forest cover maps for the reserve. Because the output of VCT disturbance year map contains 7 classes, to produce forest cover product, we needed to aggregate the 7 classes into forest and non-forest, two classes.

Table 2. Definition and aggregation of the forest disturbance map.

Code	Class Description in VCT Model	Aggregated Class
0	Background	Abandoned
1	Persisting non-forest	Non-forest
2	Persisting forest	Forest
3	Persisting water	Non-forest
4	Probable forest with recent disturbance	Forest
5	Disturbed in this year	Non-forest
6	Post-disturbance non-forest	Non-forest

shows the detailed criteria for the aggregation.

2.4. Development and Validation of NDVI_DR Model

To support subsequent AGB modelling by forest type, we must create a time-series of forest type distribution data for the reserve. The growth of vegetation is affected by various factors, such as temperature and precipitation, and presents different growth characteristics in different periods. Through analyzing the difference in vegetation spectral characteristics in different growth periods, different vegetation types can be effectively distinguished [40,41]. This work used the remote sensing images collected both in growing season (summer) and winter season of the same year to construct a new NDVI-based index to distinguish coniferous forest from broad-leaved forest. Equation (1) illustrates the specifics of the new index, NDVI-DR. Specifically, based on the 2011 forest resource inventory and planning data, we first randomly selected 200 coniferous forest stands and 200 broad-leaved forest stands, respectively, and, based on the gravity center of each stand, its 8 adjacent pixels (8 neighborhoods) in 8 directions of the central pixel were jointly considered. Thus, the average NDVI value of the 9 pixels (central pixel plus 8 adjacent pixels) was extracted as the modified NDVI value of the central pixel to minimize the potential locational errors between summer image and corresponding winter image. By using the NDVI_S (NDVI value in summer) and NDVI_W (NDVI value in winter) of the same location, a new image index,$\mathrm{NDVI}\_\mathrm{DR}$, was calculated by following Equation (1). The change patterns of NDVI_S, NDVI_W and NDVI_DR for different central pixels were plotted in Figure 2. Figure 2 shows that the NDVI values of coniferous forest pixels and broad-leaved forest pixels in summer (NDVI_S) are higher than those in winter (NDVI_W). The mean NDVI of coniferous forest is 0.742 in summer and 0.617 in winter. The mean NDVI of broad-leaved forest is 0.571 in summer and 0.209 in winter. The NDVI values of coniferous forest in summer have little differences compared with those in winter, while the NDVI values of broad-leaved forest in summer have greater differences compared with those in winter. According to the above observations, a new NDVI_DR vegetation index was constructed following Equation (1). In this equation, the denominator is the difference between NDVI in summer and NDVI in winter, and the numerator is the NDVI value in winter.

$$\mathrm{NDVI}\_\mathrm{DR}=\left(\mathrm{NDVI}\_\mathrm{S}-\mathrm{NDVI}\_\mathrm{W}\right)/\left(\mathrm{NDVI}\_\mathrm{W}\right)$$

(1)

NDVI_DR can reflect the variation intensity of NDVI value in summer and winter. Figure 2 shows that, in coniferous forest, the value of NDVI_DR is generally less than or equal to 0.4, while, in broad-leaved forest, the value of NDVI_DR is greater than or equal to 0.5. Thus, the thresholds of 0.4 and 0.5 for NDVI_DR were determined as the final classification criteria for forest type identification.

To validate the effectiveness of NDVI_DR thresholding model, we used the 2013 Landsat 8 OLI images coupled with the 2012 vegetation distribution map and adopted fully the same process to extract the values of NDVI_DR for those 400 locations or pixels. After visually interpreting the 2012 vegetation distribution map, the forest types of those 400 pixels remaining unchanged were doubly confirmed. The extracted change patterns of those 400 locations based on the 2013 Landsat OLI image were summarized in Figure 3. Obviously, the proposed NDVI_DR thresholding model remains stable; the thresholds of 0.4 and 0.5 are still effective in separating coniferous forest from broad-leaved forest. Thus, the following rules in Equation (2) were used to distinguish coniferous forest from broad-leaved forest by using Landsat observations in the current work. Equation (2) was written as:

$$\{\begin{array}{c}\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\_\mathrm{D}\mathrm{R}\le 0.4\to \mathrm{C}\mathrm{o}\mathrm{n}\mathrm{i}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{s} \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\hfill \\ \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\_\mathrm{D}\mathrm{R}\ge 0.5\to \mathrm{b}\mathrm{r}\mathrm{o}\mathrm{a}\mathrm{d}-\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{d} \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\hfill \end{array}$$

(2)

2.5. Forest AGB Modeling

2.5.1. Independent Variable

Based on the 2011 forest resources inventory and planning data, the per area AGB was derived as the independent variable for modelling. Specifically, according to the recorded attributes in the inventory data, including the dominant tree species, average diameter at breast height (DBH), average height, area and number of stems of each stand, the single tree-level AGB was calculated by using the allometric growth equations by tree species (

Table 3. Biomass allometric growth equations for the major tree species in Anhui Province.

Tree Species	Aboveground Biomass Formula
Cedarwood	${\mathrm{W}}_{\mathrm{T}}{=\mathrm{W}}_{\mathrm{S}}{+\mathrm{W}}_{\mathrm{B}}{+\mathrm{W}}_{\mathrm{L}}=0{.00849(\mathrm{D}}^2{\mathrm{H})}^{1.107230}+0{.00175(\mathrm{D}}^2{\mathrm{H})}^{1.091916}+0{.00071\mathrm{D}}^{3.88664}$
Oak	${\mathrm{W}}_{\mathrm{T}}{=\mathrm{W}}_{\mathrm{S}}{+\mathrm{W}}_{\mathrm{B}}{+\mathrm{W}}_{\mathrm{L}}=0{.00888(\mathrm{D}}^2{\mathrm{H})}^{1.08}+0{.01(\mathrm{D}}^2{\mathrm{H})}^{0.90}+0{.00378(\mathrm{D}}^2{\mathrm{H})}^{0.94}$
Larch	${\mathrm{W}}_{\mathrm{T}}{=\mathrm{W}}_{\mathrm{S}}{+\mathrm{W}}_{\mathrm{B}}{+\mathrm{W}}_{\mathrm{L}}=0{.099496(\mathrm{D}}^2{\mathrm{H})}^{0.786530}+0{.098620(\mathrm{D}}^2{\mathrm{H})}^{0.598367}+0{.294136(\mathrm{D}}^2{\mathrm{H})}^{0.357506}$
Masson pine	${\mathrm{W}}_{\mathrm{T}}=0{.01672(\mathrm{D}}^2{\mathrm{H})}^{0.8559}$
Sclerophyllous broad-leaved forest	${\mathrm{W}}_{\mathrm{T} }=0{.07112(\mathrm{D}}^2{\mathrm{H})}^{0.910359078}$
Soft-leaved broad-leaved forest	${\mathrm{W}}_{\mathrm{T}}{=\mathrm{W}}_{\mathrm{S}}{+\mathrm{W}}_{\mathrm{B}}{+\mathrm{W}}_{\mathrm{L}}=0{.012541(\mathrm{D}}^2{\mathrm{H})}^{1.144}+0{.004786(\mathrm{D}}^2{\mathrm{H})}^{1.006}+0{.047180(\mathrm{D}}^2{\mathrm{H})}^{0.769}$

Note: ${\mathrm{W}}_{\mathrm{T}}$ is the total AGB; ${\mathrm{W}}_{\mathrm{S}}$ is the trunk biomass; ${\mathrm{W}}_{\mathrm{B}}$ is the branch biomass; ${\mathrm{W}}_{\mathrm{L}}$ is the leaf biomass; $\mathrm{D}$ is the DBH of trees; $\mathrm{H}$ is the height of trees. The AGB formula is derived from the main technical provisions of China’s National Forest Inventory.

), and then the total stand-level AGB was derived by summing up the AGB of each tree in the stand, followed by the calculation of per unit area AGB (t/ha) via dividing the total stand-level AGB by the stand area. To match the pixel size of Landsat image in the AGB modelling, it was necessary to convert the per unit area AGB into the pixel-level AGB, with the unit of t/900 m². After this conversion, all the field AGB measurements were used as the independent variable for model training (80%) and validation (20%) purposes to facilitate the modelling.

2.5.2. Development of Dependent Variables

In this work, four types of modelling features, including the original band transformations, vegetation indices, textural measures and terrain variables, were extracted as the potential dependent variables to support the establishment of AGB inversion model.

The Landsat multi-spectral imagery has abundant spectral information and different bands have different levels of correlation to AGB. In order to remove the redundant information among bands and to screen out those comprehensive features highly relating to AGB, this analysis adopted the KT transform to generate three new orthogonal features with explicit physical implications, named TCB (Brightness), TCG (Greenness), TCW (Wetness) [42]. Meanwhile, the TCD (Distance) [43] and TCA (Angle) [44] indices that reflect vegetation coverage and tree growth status [45] were also developed based on the three features (Table 4).

Vegetation index can enhance vegetation signature and accurately reflect vegetation growth and distribution [46]. The ratio vegetation index (RVI) [47] can eliminate the influence of terrain and shadow on vegetation analysis, and the normalized difference vegetation index (NDVI) can minimize the effects of atmosphere on vegetation and characterize vegetation density and vigor, and both indices show a good correlation with AGB [48,49]. Additionally, the NDVI_C index, making use of its better resistance of short-wave infrared band to atmospheric condition changes, can unify different coverage types and reduce the influence of forest background signals [50,51]. These indices were derived by following the formula in Table 4.

Table 4. Vegetation indices used in the analysis and their calculation formulas.

Index	Formula
NDVI [48]	$\mathrm{NDVI}=\frac{{\mathsf{\rho }}_{\mathrm{nir}}-{\mathsf{\rho }}_{\mathrm{red}}}{{\mathsf{\rho }}_{\mathrm{nir}}+{\mathsf{\rho }}_{\mathrm{red}}}$
NDVIC [50]	${\mathrm{NDVI}}_{\mathrm{C}}=\frac{{\mathsf{\rho }}_{\mathrm{nir}}-{\mathsf{\rho }}_{\mathrm{red}}}{{\mathsf{\rho }}_{\mathrm{nir}}+{\mathsf{\rho }}_{\mathrm{red}}}\times \left(1-\frac{{\mathsf{\rho }}_{\mathrm{swir}}-{\mathsf{\rho }}_{{\mathrm{swir}}_{\min }}}{{\mathsf{\rho }}_{{\mathrm{swir}}_{\max }}-{\mathsf{\rho }}_{{\mathrm{swir}}_{\min }}}\right)$
RVI43 [47]	${\mathrm{RVI}}_{43}=\frac{{\mathsf{\rho }}_{\mathrm{nir}}}{{\mathsf{\rho }}_{\mathrm{red}}}$
RVI54 [47]	${\mathrm{RVI}}_{54}=\frac{{\mathsf{\rho }}_{\mathrm{swir}}}{{\mathsf{\rho }}_{\mathrm{nir}}}$
NDMI [52]	$\mathrm{NDMI}=\frac{{\mathsf{\rho }}_{\mathrm{nir}}-{\mathsf{\rho }}_{\mathrm{swir}}}{{\mathsf{\rho }}_{\mathrm{nir}}+{\mathsf{\rho }}_{\mathrm{swir}}}$
mNDWI [52]	$\mathrm{mNDWI}=\frac{{\mathsf{\rho }}_{\mathrm{green}}-{\mathsf{\rho }}_{\mathrm{swir}}}{{\mathsf{\rho }}_{\mathrm{green}}+{\mathsf{\rho }}_{\mathrm{swir}}}$
TCD [43]	$\mathrm{TCD}=\sqrt{{\mathrm{TCB}}^2+{\mathrm{TCG}}^2}$
TCA [44]	$\mathrm{TCA}=\arctan \left(\frac{\mathrm{TCG}}{\mathrm{TCB}}\right)$

Note: ${\mathsf{\rho }}_{\mathrm{nir}}$ is the spectral reflectance of near infrared band. ${\mathsf{\rho }}_{\mathrm{red}} \mathrm{is} \mathrm{the} \mathrm{spectral} \mathrm{reflectance} \mathrm{of} \mathrm{red} \mathrm{band}.$ ${\mathsf{\rho }}_{\mathrm{swir}}$ is the spectral reflectance of short-wave infrared band. ${\mathsf{\rho }}_{\mathrm{green}}$ is the spectral reflectance of green band.

Table 4. Vegetation indices used in the analysis and their calculation formulas.

Index	Formula
NDVI [48]	$\mathrm{NDVI}=\frac{{\mathsf{\rho }}_{\mathrm{nir}}-{\mathsf{\rho }}_{\mathrm{red}}}{{\mathsf{\rho }}_{\mathrm{nir}}+{\mathsf{\rho }}_{\mathrm{red}}}$
NDVIC [50]	${\mathrm{NDVI}}_{\mathrm{C}}=\frac{{\mathsf{\rho }}_{\mathrm{nir}}-{\mathsf{\rho }}_{\mathrm{red}}}{{\mathsf{\rho }}_{\mathrm{nir}}+{\mathsf{\rho }}_{\mathrm{red}}}\times \left(1-\frac{{\mathsf{\rho }}_{\mathrm{swir}}-{\mathsf{\rho }}_{{\mathrm{swir}}_{\min }}}{{\mathsf{\rho }}_{{\mathrm{swir}}_{\max }}-{\mathsf{\rho }}_{{\mathrm{swir}}_{\min }}}\right)$
RVI43 [47]	${\mathrm{RVI}}_{43}=\frac{{\mathsf{\rho }}_{\mathrm{nir}}}{{\mathsf{\rho }}_{\mathrm{red}}}$
RVI54 [47]	${\mathrm{RVI}}_{54}=\frac{{\mathsf{\rho }}_{\mathrm{swir}}}{{\mathsf{\rho }}_{\mathrm{nir}}}$
NDMI [52]	$\mathrm{NDMI}=\frac{{\mathsf{\rho }}_{\mathrm{nir}}-{\mathsf{\rho }}_{\mathrm{swir}}}{{\mathsf{\rho }}_{\mathrm{nir}}+{\mathsf{\rho }}_{\mathrm{swir}}}$
mNDWI [52]	$\mathrm{mNDWI}=\frac{{\mathsf{\rho }}_{\mathrm{green}}-{\mathsf{\rho }}_{\mathrm{swir}}}{{\mathsf{\rho }}_{\mathrm{green}}+{\mathsf{\rho }}_{\mathrm{swir}}}$
TCD [43]	$\mathrm{TCD}=\sqrt{{\mathrm{TCB}}^2+{\mathrm{TCG}}^2}$
TCA [44]	$\mathrm{TCA}=\arctan \left(\frac{\mathrm{TCG}}{\mathrm{TCB}}\right)$

Note: ${\mathsf{\rho }}_{\mathrm{nir}}$ is the spectral reflectance of near infrared band. ${\mathsf{\rho }}_{\mathrm{red}} \mathrm{is} \mathrm{the} \mathrm{spectral} \mathrm{reflectance} \mathrm{of} \mathrm{red} \mathrm{band}.$ ${\mathsf{\rho }}_{\mathrm{swir}}$ is the spectral reflectance of short-wave infrared band. ${\mathsf{\rho }}_{\mathrm{green}}$ is the spectral reflectance of green band.

The texture of remote sensing image refers to the recurring primitives or elements and their arrangement rules in the image, which is a unity of local variability and spatial correlation [53,54,55]. In this paper, 6 multispectral bands, NDVI_C and RVI₅₄ were used as the inputs for texture calculation by using the gray level co-occurrence matrix (GLCM) method [56]. Specifically, a 3 × 3 window size was selected, and the moving directions at 0°, 45°, 90° and 180° were considered, respectively, and the average of the four directions was taken as the final texture analysis result. To compensate for the potential limitations of GLCM method, the wavelet multi-scale decomposition was also implemented by programming in MATLAB environment to extract the high-frequency vertical, horizontal and diagonal details images as new textural features [57]. Firstly, the first principal component of the 2011 Landsat multi-spectral images was decomposed by using a biorthogonal wavelet base function in a three-layer recursive manner. Thus, another 9 images of details were obtained as new textural features to support the modeling. The GLCM-based textures and their calculation formula were summarized in

Table 5. Formulas for calculating texture index based on GLCM.

Texture Index	Formula
Mean, (ME)	$\mathrm{ME}=\sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{i} \times \mathrm{P}}_{\mathrm{i},\mathrm{j}}$
Variance, (VA)	$\mathrm{VA}=\sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{i} \times \mathrm{P}}_{\mathrm{i},\mathrm{j}}{\left(\mathrm{i} - \mathrm{ME}\right)}^2$
Homogeneity, (HO)	$\mathrm{HO}=\sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}\mathrm{i} \times \frac{{\mathrm{P}}_{\mathrm{i},\mathrm{j}}}{1+{\left(\mathrm{i}-\mathrm{j}\right)}^2}$
Contrast, (CO)	$\mathrm{CO}=\sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{i} \times \mathrm{P}}_{\mathrm{i},\mathrm{j}}{\left(\mathrm{i} - \mathrm{j}\right)}^2$
Dissimilarity, (DI)	$\mathrm{DI}=\sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{i} \times \mathrm{P}}_{\mathrm{i},\mathrm{j}}\left|\mathrm{i} - \mathrm{j}\right|$
Entropy, (EN)	$\mathrm{EN}=\sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{i} \times \mathrm{P}}_{\mathrm{i},\mathrm{j}}\left({-\mathrm{lnP}}_{\mathrm{i},\mathrm{j}}\right)$
Second Moment, (SM)	$\mathrm{SM}=\sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{i} \times \mathrm{P}}_{\mathrm{i},\mathrm{j}}{}^2$
Correlation (CR)	$\mathrm{CR}=\sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{i} \times \mathrm{P}}_{\mathrm{i},\mathrm{j}}\left[\frac{\left(\mathrm{i} - \mathrm{ME}\right)\left(\mathrm{j} - \mathrm{ME}\right)}{\sqrt{{\mathrm{VA}}_{\mathrm{i}}{•\mathrm{VA}}_{\mathrm{j}}}}\right]$

Note: ${\mathrm{P}}_{\mathrm{i},\mathrm{j}}=\frac{{\mathrm{V}}_{\mathrm{i},\mathrm{j}}}{{\sum }_{\mathrm{i},\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{V}}_{\mathrm{i},\mathrm{j}}};{ \mathrm{V}}_{\mathrm{i},\mathrm{j}}$ represents the pixel brightness value at the position of row i and column j; N is the size or dimension of the moving window when the texture index is calculated.

[58].

Terrain features influence the distribution patterns of AGB to some extent [59]. Therefore, elevation, slope and aspect extracted from the 12.5 m resolution DEM were also considered as the modeling factors. These terrain factors were resampled to 30 m resolution by implementing the bi-linear interpolation resampling technique to match the Landsat pixel size.

2.5.3. Correlation Analysis

In this paper, 6 multispectral bands, 3 topographic factors, 6 vegetation indices, 5 KT transform and 73 texture features were selected as potential variables for biomass prediction. In order to ensure the accuracy of the model and reduce the workload, correlation analysis was needed to select the best combination of variables [60]. Assuming that there is no correlation or weak correlation between the selected variables and a linear relationship with the dependent variable, AGB, we could automatically select the variables by regression analysis, and this method was usually used in previous studies [61,62]. This analysis first used random forest importance ranking to select those important variables and then eliminated variables with high correlation between variables through Pearson correlation analysis [63]. The random forest package was run in R environment, and the characteristic factors with high correlation to biomass of broad-leaved forest and coniferous forest were, respectively, screened out as modeling variables.

2.5.4. SGB-Based AGB Modelling and Its Time Extrapolation

Stochastic gradient boosting (SGB) algorithm proposed by Friedman was adopted to estimate AGB in this work [64]. SGB takes into account the advantages of both boosting algorithm and bagging algorithm and has been widely used in regression and classification tasks [65]. It can avoid the problem of long calculation time due to large amount of data and can also improve the prediction accuracy, with better robustness to overfitting [33]. The “learning rate” parameter in the SGB algorithm determines the growth rate of modeling complexity. Generally, a smaller learning rate means that more regression trees will be generated and the contribution of each tree to the whole forest will be weaker and the modeling performance will be better [66]. The “depth of interaction” parameter determines the splitting number of each tree. This parameter value represents the number of nodes in each decision tree, and the maximum value is usually set to 10. Other important parameters include: tree complexity, shrinkage, distribution function, training ratio, etc. The modeling process was implemented by using the “caret” package in R environment. Through multiple comparisons, it was found that the optimal combination of parameters in SGB modeling was set as: interaction depth at 3, shrinkage at 0.01, ntrees at 500 and n.minobsinnode at 9. Thus, the 2011 AGB map of the reserve was produced from implementing the above-identified SGB parameters by forest types, followed by a spatial overlay analysis of both the estimated coniferous forest AGB and the predicted broad-leaved forest AGB.

At present, there are two main methods to observe long-term forest AGB changes using remote sensing images. The first is to build separate varying relationships between remote sensing images and corresponding forest AGB field sample measurements in different years. This method is relatively accurate, but it is subject to the constraints of realistic conditions, such as lack of historical field AGB measurement data. The second method assumes that there is a relatively stable relationship between forest AGB and remote sensing images over time. Through the relative radiometric normalization operations between the base image (e.g., the 2011 Landsat TM image of the current work) and the target images in different years, followed by the development of a fully same set of features as the base image based on the normalized target images, the relationship constructed for the base year is extended to other target years to retrieve varying AGB patterns in the same study area. The second method has been more widely used in estimating forest AGB thanks to its lower cost and stronger operability in a long-term monitoring period [67,68]. Here, the 2011 Landsat TM image was set as the base image, and Landsat TM/ETM+/OLI images acquired in 1987, 1992, 1997, 2002, 2007, 2017 and 2020, as the target images, were radiometrically normalized, respectively, band by band using the automated weighted invariant points (WIP) method [69]. Based on these normalized target images, a fully same set of features as the base year was developed. Then, the relationship created for the base year was extended to other target years to create a time-series of AGB maps for the reserve.

2.6. Validation Method

2.6.1. Forest Distribution Verification

According to the historical documentations of the reserve, forest was the largest land cover type in the reserve, occupying an area proportion of about 90%. Thus, to validate the accuracy of forest distribution products mapped from VCT, the stratified random sampling method was implemented. First, based on the 2011 forest resources inventory and planning data and the 2012 vegetation distribution map, 900 points were randomly generated in the forest area and 100 points were randomly generated in the non-forest area for the classifications in other years. Then, visually interpreting the corresponding year’s Google Earth images was conducted to gain the ground truths to validate the classifications. For the year of 1997 or earlier, no Google Earth high resolution images were available, so the original Landsat images were directly visually interpreted. Finally, by comparing the classifications and the interpreted results, the overall accuracy and kappa coefficient were derived to evaluate the classification accuracy. The kappa coefficient calculation formula was as follows:

$$\mathrm{kappa}=\frac{{\mathrm{P}}_{\mathrm{A}}-{\mathrm{P}}_{\mathrm{e}}}{1-{\mathrm{P}}_{\mathrm{e}}}$$

(3)

where ${\mathrm{P}}_{\mathrm{A}}$ refers to the relative observed agreement among raters, and ${\mathrm{P}}_{\mathrm{e}}$ is the hypothetical probability of chance agreement.

2.6.2. Forest Type Distribution Verification

The broad-leaved forest dominated the reserve historically. Similarly, the stratified random sampling method was used to verify the forest type classification accuracy. Based on the 2011 forest resources inventory and planning data and the 2012 vegetation distribution map, 700 points were randomly generated in the broad-leaved forest area and 300 points were randomly generated in the coniferous forest area. Then, visually interpreting the corresponding year’s Google Earth images was conducted to gain the ground truths to validate the classifications.

2.6.3. Forest AGB Modeling Verification

Twenty percent of the pixel-level AGB field measurements, as an independent dataset, were used to evaluate the SGB-based modelling of AGB by calculating the validation determination coefficient R² and root mean square error (RMSE). R² indicates the variance explanation degree of dependent variables to independent variable in the model. Generally speaking, larger R² and smaller RMSE mean higher prediction accuracy that the model has. The calculation formulas of R² and RMSE were shown in Equations (4) and (5):

$${\mathrm{R}}^2=1-\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\widehat{{\mathrm{y}}_{\mathrm{i}}}\right)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\stackrel{-}{\mathrm{y}}\right)}^2}$$

(4)

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\widehat{{\mathrm{y}}_{\mathrm{i}}}\right)}^2}{\mathrm{n}}}$$

(5)

where $\widehat{{\mathrm{y}}_{\mathrm{i}}}$ indicates the model-predicted AGB, ${\mathrm{y}}_{\mathrm{i}}$ represents the measured AGB, $\stackrel{-}{\mathrm{y}}$ indicates the average value of the measured AGB and $\mathrm{n}$ means the sample size.

Based on the 2011 forest resources inventory data in 2011, 200 broad-leaved forest stands and 200 coniferous forest stands were randomly selected, of which 120 stands were used for modeling and 80 for independent validation. Additionally, to highlight the potential superior performance of separate modelling by forest types, these 400 stands were combined to construct a mixed AGB model without differentiating forest types, of which 240 stands were used for modeling and 160 for validation. Because the historical field records of AGB in the reserve were unavailable, we were unable to validate the AGB predictions produced from extrapolating the 2011 AGB models to other years.

3. Results

3.1. Forest Distribution Mapping and Validation

Table 6. Confusion matrix of forest extraction accuracy verification.

	Predicter Results	Forest	Non-Forest	Total
Actual Results
Forest		863	16	879
Non-Forest		37	84	121
Total		900	100	1000
OA = 0.947			Kappa = 0.731

shows the confusion matrix of the 2020 forest classifications. The results showed that an overall accuracy of 94.7%, accompanied by a kappa coefficient of 0.73, were observed in the 2020 forest classifications.

Table 7. Verification results of forest extraction accuracy.

Year	Overall Accuracy (%)	Kappa Coefficient
1987	93.2%	0.66
1992	93.1%	0.66
1997	92.8%	0.64
2002	90.5%	0.56
2007	93.9%	0.70
2011	93.1%	0.67
2017	90.2%	0.55
2020	94.7%	0.73

exhibits the verification results in other years. As a whole, the overall accuracy of the VCT-based forest classification products was over 90% and the kappa coefficient above 0.6, which meant that VCT produced reliable time-series forest distribution maps.

Figure 4 shows the spatio-temporal distribution patterns of the forest in the Yaoluoping National Nature Reserve during the period of 1987 to 2020. Obviously, the forest absolutely dominated the reserve in all the years and its areal proportion dropped from 92.92% in 1987 to 87.15% in 1997, then gradually fluctuated to 90.94% in 2020 (Figure 5). The lowest forest area was observed in 1997, and then forest area showed an upward fluctuating trend. The non-forest type only occupied a small area, principally distributed in the northern and central portions of the reserve (Figure 4), along with the valleys and the gentle-slope lands at relatively low altitudes in the study area.

3.2. Forest Type Distribution Mapping

Table 8. Verification results of forest type classifications derived from NDVI_DR thresholding.

Year	Overall Accuracy (%)	Kappa Coefficient
1987	94.5%	0.87
1992	94.0%	0.85
1997	94.2%	0.86
2002	92.5%	0.82
2007	93.5%	0.84
2011	94.0%	0.85
2017	92.0%	0.81
2020	94.8%	0.87

displays the validation accuracy statistics of the forest type classifications derived from the NDVI_DR thresholding. The overall accuracies at 94.5%, 94.0%, 94.2%, 92.5%, 93.5%, 94.0%, 92.0% and 94.8%, with corresponding kappa coefficient values at 0.87, 0.85, 0.86, 0.82, 0.84, 0.85, 0.81 and 0.87, were observed during the period of 1987 to 2020, respectively. Among them, the imaging times of the 2002 and 2017 summer images were on September 28 and October 7, respectively, which fell into late fall season and were not in the peak season of plant growth. This would lead to actual changes in the vegetation spectral reflectance, accordingly resulting in a relatively low accuracy of forest type classifications in these two years. However, overall, the DNVI_DR thresholding model provided a relatively high accuracy of forest type classifications in the reserve (Table 8).

Figure 6 shows the spatio-temporal variations in the forest type distribution. Generally, broad-leaved forest dominated the reserve in space, with an apparent large-scale continuous distribution in all the years (Figure 6). Its areal share remained relatively stable (above 80% in all the years), with a small variation (Figure 7). Coniferous forest had a dispersing distribution over the reserve, mainly distributed in the northern and south-eastern portions of the reserve (Figure 6), and its areal share gradually decreased from 8.5% in 1987 to the lowest value at about 4.7% in 2011 (Figure 7), and then the proportion of coniferous forest area increased steadily. By 2020, the proportion slightly exceeded that in 1987.

3.3. Biomass Estimation Results

3.3.1. Variable Selection Results

Table 9. Characteristic variables with significant correlation to broad-leaved forest AGB.

Characteristic	Correlation Index	Characteristic	Correlation Index
B1	−0.152 *	RVI54	0.294 **
B2	−0.314 **	NDMI	0.287 *
B3	−0.221 *	mNDWI	0.267 *
B4	−0.265 *	B5_mean	−0.473 **
B5	−0.469 **	B7_mean	−0.473 **
B7	−0.486 **	RVI54_mean	0.211 *
TCW	−0.478 **	NDVIC _mean	0.452 **
TCB	−0.406 **	B2_ Correlation	0.166 *
TCG	−0.263 **	RVI54_ Correlation	0.154 *
TCD	−0.341 **	DBH	0.326 **
NDVI	0.305 **	bio2.8	0.196 *
NDVIC	0.466 **

Note: B1, B2, B3, B4, B5, B7 mean the blue, green, red, NIR, SWIR 1, SWIR 2 band of Landsat 5 and Landsat 8 separately. * means significant at the level of 0.05 (the confidence level is 95%); ** means significant at the level of 0.01 (the confidence level is 99%).

,

Table 10. Characteristic variables with significant correlation to coniferous forest AGB.

Characteristic	Correlation Index	Characteristic	Correlation Index
B1	−0.138 *	RVI54	0.324 **
B2	−0.301 **	NDMI	0.246 *
B3	−0.191 *	mNDWI	0.207 **
B4	−0.224 *	NDVI	0.312 **
B5	−0.435 **	NDVIC	0.459 **
B7	−0.477 **	bio2.8	0.186 *
TCD	−0.264 **

Note: * means significant at the level of 0.05; ** means significant at the level of 0.01.

and

Table 11. Characteristic variables with significant correlation to combined AGB of both forest types.

Characteristic	Correlation Index	Characteristic	Correlation Index
B2	−0.173 *	RVI54	0.163 *
B4	−0.107 *	NDMI	0.142 *
B5	−0.227 *	mNDWI	0.101 *
B7	−0.248 **	NDVI	0.152 *
TCD	−0.134 *	NDVIC	0.279 **

Note: * means significant at the level of 0.05; ** means significant at the level of 0.01.

show the combinations of variables determined by using random forest importance ranking coupled with Pearson correlation analysis, most suitable for broad-leaved forest, coniferous forest and the combination of both in AGB modeling. As shown in Table 9, the average texture information of SWIR and NDVI_C had a high correlation with the broad-leaved forest AGB, while the absolute value of the correlation coefficient of other texture information was less than 0.3. As can be seen from Table 10, short-wave infrared, RVI₅₄, NDVI and NDVI_C were highly correlated with coniferous forest AGB, and the correlation between texture features and coniferous forest AGB was very low; thus, they were unable to be selected. As can be seen from Table 11, when considering the combined AGB of the coniferous forest and broad-leaved forest together, the number of selected modeling variables was substantially less than that when distinguishing between the forest types, and the correlation coefficients were also apparently lower. Finally, according to the ranking of the correlation coefficient values, 10 modeling factors were selected for broad-leaved forest, coniferous forest and combined forest types together, respectively, to consider the operational practicability. The factors selected for broad-leaved forest were B7, B5, TCW, TCB, TCD and B5_mean, B7_mean, NDVI_C _mean and NDVI_C; RVI₅₄, NDVI_C, B2, B5, B4, B7, TCD, NDVI, mNDWI and NDMI for coniferous forest; B2, B4, B5, B7, TCD, RVI₅₄, NDMI, mNDWI, NDVI and NDVI_C for combined AGB modelling of both forest types.

3.3.2. Modeling Accuracy Evaluation

Figure 8, Figure 9 and Figure 10 show the modeling and validation accuracy of broad-leaved forest, coniferous forest and combined forest types using SGB, respectively. The modeling and validation R² of broad-leaved forest AGB were at 0.68 and 0.63, and the corresponding RMSEs were at 7.53 and 11.19 t/hm², respectively. The modeling and validation R² of coniferous forest biomass were at 0.71 and 0.61, and the corresponding RMSEs were 4.46 and 14.27 t/hm², respectively. The modeling and validation R² of combined biomass were at 0.54 and 0.51, with the corresponding RMSEs at 18.91 and 20.47 t/hm², respectively. Obviously, modelling AGB by forest type was more accurate than modelling AGB without differentiating forest types.

3.3.3. AGB Extrapolation of SGB Model

In this work, the 2011 well-established and tested SGB-based AGB models for broad-leaved forest and coniferous forest were extrapolated to create a time-series of AGB based on the forest type maps created in Section 3.2 during the period of 1987 to 2020 (Figure 11). The total biomass and mean values of different forest types in each year were calculated, as shown in

Table 12. Statistics of the mapped AGB in each year.

Year	Forest		Coniferous Forest		Broad-Leaved Forest
	Mean (t/hm2)	Summation (10,000 Tons)	Mean (t/hm2)	Summation (10,000 Tons)	Mean (t/hm2)	Summation (10,000 Tons)
1987	57.37	70.57	62.27	6.77	56.68	60.20
1992	61.56	75.72	64.84	5.41	60.12	65.20
1997	64.38	79.19	67.53	5.59	63.74	64.72
2002	70.14	86.27	72.03	7.37	69.45	69.76
2007	73.21	90.05	74.71	5.45	72.37	76.23
2011	75.52	92.89	76.43	4.08	74.72	80.84
2017	76.23	93.76	78.98	6.97	75.46	77.59
2020	78.85	96.99	80.45	8.88	78.85	81.64

. In terms of temporal distribution, the forest AGB in the reserve showed a gradual increasing trend in the past 30 years. The average AGB values in the eight years were at 57.37, 61.56, 64.38, 70.14, 73.21, 75.52, 76.23 and 78.85 t/hm², respectively, which showed that the AGB of Yaoluoping National Nature Reserve was increasing continuously, and the growth rate was the fastest during the period of 1992 to 2002. As for different forest types, although the total AGB of coniferous forest was much lower than that of broad-leaved forest, the average biomass of coniferous forest was higher than that of broad-leaved forest (Table 12).

The spatial distribution of the AGB in Figure 11 shows that the areas with low AGB values were mainly located near agricultural lands, buildings and water bodies in river valleys, and took on apparent linear features in the northern and central portions of the reserve, whereas high AGB values mainly occurred in the western and north-eastern areas of the reserve (Figure 11).

4. Discussion

4.1. NDVI_DR Thresholding Model

Accurate identification of forest types is of great significance to forestry development planning and forestry policy formulation, helps to reliably evaluate the successive stage and trend of a specific forest ecosystem [3,4] and potentially improves the modelling accuracy of remote-sensing-based AGB estimation due to being able to separately model coniferous forest and broad-leaved forest [14,15,16]. In this paper, the NDVI_DR vegetation index that makes use of phenological or seasonal differences in the spectral signature of evergreen forest and deciduous forest was developed to classify forest types in the study area based on long-term time-series Landsat images, with an overall accuracy of above 92% and a kappa coefficient of about 0.85. The classification accuracy of the current NDVI_DR thresholding model is higher than other similar studies using Landsat for forest type classification. For example, Li et al. [70] used three machine learning approaches, including decision trees, random forest and support vector machines and Landsat images, to classify local forest communities at the Huntington Wildlife Forest (HWF). Among them, the SVM had the highest accuracy, with an overall accuracy of 88.2% and kappa coefficient of 0.793. Hill et al. [71] used two methods of low-pass spatial filtering to reduce the local spectral variation and image segmentation to implement supervised classifications of forest types in Peruvian Amazonia from Landsat TM data and gained an overall accuracy of about 90%. In their studies, they constructed various vegetation indices to express the spectral characteristics differences in different forest types, and they found that the spectral differences in different forest types characterized by diverse vegetation indices were not adequately separable. In contrast to these existing investigations, our NDVI_DR thresholding model does not require the development of diverse image features, e.g., spectral indices, image textures or contextual information; more importantly, there is no need to accurately tune the parameters of the classification models or algorithms. Conversely, our NDVI_DR model just calculates a derivative of a forest pixel’s NDVI in the summer and corresponding NDVI in the winter to reflect the seasonal differences in the spectral signature of forest types, and it specifies stable thresholds to the derivative (less than 0.4 for coniferous forest and greater than or equal to 0.5 for broad-leaved forest, Figure 2 and Figure 3) to classify forest types. Obviously, our model is efficient and easily implemented compared to the existing methods, and it substantially improves the identification accuracy of forest types by considering the seasonal differences of different forest types, which is in agreement with Zhu’s [72] research results. Dong et al. [73] found that using seasonal time-series data has the potential to improve the accuracy for monitoring forest attributes, whereas the separability and stability of the NDVI_DR thresholding model have just been tested in the subtropical forest ecosystem by Landsat TM and OLI sensors, so its robustness in other ecoregions and other similar sensors, such as Sentinel-2 MSI, should be continuously verified in the near future to doubly confirm its generalization or popularization.

4.2. Forest AGB Modeling

At present, the main methods used in forest biomass prediction include random forest, support vector machine and a multiple linear regression model [20,21]. Nguyen et al. [74] developed a random-forest (RF)-based kNN model to produce annual maps of AGB from 1988 to 2017 for over 7.2 million ha of forests in Victoria, Australia. The modeling R² is between 0.37 and 0.59, and the RMSE is between 104.7 and 168.5 t/hm². Main-Knorn et al. [75] obtained the AGB from 1985 to 2010 through RF models based on Landsat time-series images and field data, showing an RMSE of 41.3 t/hm². Compared with the previous methods, the SGB model adopted in this paper has better robustness for outliers, inaccurate data, missing values and unbalanced datasets, and has relatively stable estimation results [64]. Dube [33] found that, when Landsat series images were selected as research data, the SGB algorithm had higher accuracy in forest biomass estimation than the random forest algorithm regardless of texture features, spectral features or if both of them were used. In addition, most existing forest AGB modelling and mapping studies did not distinguish between forest types but considered forests as a whole to retrieve the forest aboveground biomass in the entire region [76,77]. However, Fassnacht [78] actually found that there were fundamental differences in the NIR reflectance between coniferous forests and broad-leaved forests, and hardwood canopies could reflect 50% more in NIR than pine canopies due to different cellulose compactness or structures in their leaves. Thus, separate AGB modeling of different forest types may more adequately capture the respective variances in canopies’ signatures of coniferous forest and broad-leaved forest. Actually, Figure 8, Figure 9 and Figure 10 show the modeling R² of broad-leaved forest, coniferous forest and combined forest types at 0.68, 7.53 and 0.54, respectively, with the corresponding RMSEs at 7.53, 4.46 and 14.27 t/hm². This shows that AGB modeling by different forest types can achieve higher accuracy of AGB estimation than modelling without differentiating between forest types [14,15,16], which is consistent with Zheng’s study [79].

However, due to the unavailability of historical records of field sample plots, we were unable to reliably validate the long-term time-series AGB maps generated by extrapolating the 2011 AGB prediction models into other years in the current analysis (Figure 11), although these AGB maps can provide basic data support for evaluating the ecosystem dynamic trends of the reserve and the effectiveness of management practices. Although we have assumed that there is a relatively stable relationship between forest AGB and remote sensing images in a specific area over time, in reality, natural disasters or anthropogenic disturbances may alter the forest spatial structure, species compositional structure and age structure over time and space, which may affect this stable relationship, thus affecting the accuracy of those extrapolated AGB estimations [80,81]. In the future, more efforts should be invested to ensure a sufficient validation of the historical AGB patterns generated from satellite image archives.

4.3. Driving Factors for Forest Area and AGB Changes

Our results showed that the forest areal proportion dropped from 92.92% in 1987 to 87.15% in 1997, and from 90.11% in 2011 to 88.6% in 2017 (Figure 5). It can be seen that coniferous forest is mainly distributed in the northern and south-eastern portions of the reserve, and the increased area of the coniferous forest in the reserve is significantly more than the reduced area in the past 30 years (Figure 6). This result is in agreement with the vegetation distribution of the Yaoluoping Nature Reserve studied by Xie et al. [82].

The coniferous forest decreased mainly in the southwest, northwest and southeast of the reserve between 1987 and 2002 (Figure 6). Major forest harvesting species, such as Chinese fir and Pinus taiwanensis Hayata, grew in the southwest and southeast of the reserve. The economic sources of the local residents in Yaoluoping Reserve mainly depend on timber and agricultural production. Thus, unrestricted timber harvesting of these coniferous species led to a significant reduction in area. Bahurudeen et al. [83] found that the coniferous forest had excellent materials and high economic value and was widely used in a variety of industries, which effectively explains the rapid decline in the coniferous forest in the reserve. In addition, in 2011, a large area of coniferous forest decreased due to natural disasters, such as landslides, in the southeastern part of the reserve. The increased area of coniferous forest was mainly distributed in the valley of the Baojia River basin, and the years were mainly from 1992 to 2007. In 2001, a zero-felling quota of commercial timber was implemented, and the felling of commercial timber was basically eradicated and the overall forest felling decreased dramatically. At the same time, a large number of tree species, such as Chinese fir, have been planted and renewed in the reserve. According to the records, the natural forest of Huangliyuan forest farm has been completely updated to the existing Chinese fir forest since 1993. All the events can reasonably explain the dynamic patterns of coniferous forest in the reserve during the period of 1987 to 2020.

The biomass inversion results showed that the forest AGB in Yaoluoping National Nature Reserve increased continuously from 1987 to 2020, but the coniferous forest biomass decreased correspondingly in those disturbed years (Table 12). This indicates that human factors are the main factors affecting the forest AGB. During the periods of rapid population density and GDP growth, forests, as the main source of income and livelihood, are greatly affected, and the forest AGB decreases accordingly. Rozendaal et al. [84] showed that human activities, especially logging disturbances, had a significant impact on forest biomass. In the future development process, the local government should properly consider the carrying capacity of the forest ecosystem to population density and establish an ecological compensation mechanism combined with its own characteristic forest industries so as to realize the goal of protecting natural resources and protecting the interests of the community. Liu et al. [85] took Yaoluoping National Nature Reserve as the research object and discussed how to effectively promote the development of the national nature reserve by using ecological compensation mechanisms based on the observed problems in the development process. For the scientific protection of the reserve, we can cultivate biological resources with economic benefits and carry out activities, such as ecotourism and popularization of science, to mitigate the impacts of human demands on forest resources. Xu et al. [86] quantified the biodiversity value of Yaoluoping Nature Reserve and confirmed the ecological and economic value of the reserve, which poses a high priority of sustainable natural resources management in the reserve. In addition, some deforested areas should be classified according to the degree of deforestation and site conditions, and the vegetation restoration should be organically combined with natural regeneration and artificial intervention to shorten the time of vegetation restoration to the forest community scientifically and effectively.

4.4. Limitations and Future Improvements

Although important results were obtained from this study, the following aspects still need to be further studied: (1) Although the constructed vegetation index NDVI_DR can better distinguish between coniferous forest and broad-leaved forest, its extraction effects for mixed forest and shrubs need further verification [72]. (2) More efforts should be made to ensure a reliable validation of historical AGB maps, such as collecting measurements from those permanent sample plots possibly deployed in the reserve. (3) Although the SGB algorithm can reduce the variance and bias and has high estimation accuracy, it is very sensitive to the change in the outliers in the training samples [64]. Thus, its tuned parameters need to be optimized in later research.

5. Conclusions

In this study, a VCT model was used to generate forest cover datasets in the study area, and NDVI-DR was developed to classify the forest types. On this basis, the SGB algorithm and extrapolation were applied to retrieve the AGB in the study area from 1987 to 2020. The findings from our study can provide potential insights for long-term forest remote sensing observations, forest type classification, and accurate AGB mapping. These findings also can inform similar management agencies of a carbon accounting data basis and provide informed actions on sustainable development with high ecological interests. Based on the findings, it is concluded that:

(1)

The NDVI_DR thresholding provides an efficient and accurate classification method for distinguishing between coniferous forest and broad-leaved forest. The overall accuracy is above 92%, with a kappa coefficient above 0.8.

(2)

The 2011 forest-type-dependent stochastic-gradient-boosting-based (SGB-based) AGB estimation model achieved an independent validation R square at 0.63 and an RMSE at 11.18 t/ha for broad-leaved forest, and 0.61 and 14.26 t/ha for coniferous forest. A time-series of AGB was generated by extrapolating the 2011 AGB models to other years, and the mapped AGB showed a gradual increasing trend over the past three decades.

(3)

There is a significant correlation between human disturbance and AGB, especially irregular deforestation. Thus, we suggest that the local government should properly consider the carrying capacity of the forest ecosystem to population density and establish an ecological compensation mechanism combined with its own characteristic forest industries.