Modelling the Dynamics of Carbon Storages for Pinus densata Using Landsat Images in Shangri-La Considering Topographic Factors

Abstract

Accurate estimation of forest carbon storage is essential for understanding the dynamics of forest resources and optimizing decisions for forest resource management. In order to explore the changes in the carbon storage of Pinus densata in Shangri-La and the influence of topography on carbon storage, two dynamic models were developed based on the National Forest Inventory (NFI) and Landsat TM/OLI images with a 5-year interval change and annual average change. The three modelling methods used were partial least squares (PLSR), random forest (RF) and gradient boosting regression tree (GBRT). Various spectral and texture features of the images were calculated and filtered before modelling. The terrain niche index (TNI), which is able to reflect the combined effect of elevation and slope, was added to the dynamic model, the optimal model was selected to estimate the carbon storage, and the topographic conditions in areas of change in carbon storage were analyzed. The results showed that: (1) The dynamic model based on 5-year interval change data performs better than the dynamic model with annual average change data, and the RF model has a higher accuracy compared to the PLSR and GBRT models. (2) The addition of TNI improved the accuracy, in which R² is improved by up to 10.48% at most, RMSE is reduced by up to 7.32% at most, and MAE is reduced by up to 8.89% at most, and the RF model based on the 5-year interval change data has the highest accuracy after adding TNI, with an R² of 0.87, an RMSE of 3.82 t-C·ha⁻¹, and a MAE of 1.78 t-C·ha⁻¹. (3) The direct estimation results of the dynamic model showed that the carbon storage of Pinus densata in Shangri-La decreased in 1987–1992 and 1997–2002, and increased in 1992–1997, 2002–2007, 2007–2012, and 2012–2017. (4) The trend of increasing or decreasing carbon storage in each period is not exactly the same on the TNI gradient, according to the dominant distribution, as topographic conditions with lower elevations or gentler slopes are favorable for the accumulation of carbon storage, while the decreasing area of carbon storage is more randomly distributed topographically. This study develops a dynamic estimation model of carbon storage considering topographic factors, which provides a solution for the accurate estimation of forest carbon storage in regions with a complex topography.

1. Introduction

Forests are the mainstay of terrestrial ecosystems, which store 60% of the carbon in the terrestrial ecosystem [1]. Forests play an indispensable role in balancing and regulating CO₂ in the atmosphere with their powerful carbon sink function. Global warming is a great threat to humanity today, and the excessive emissions of greenhouse gases, of which CO₂ is an important member, is one of the main reasons for global warming. Atmospheric CO₂ concentrations have risen dramatically over the past 100 years (they have risen by approximately 110 ppm) [2], a phenomenon that has accelerated the rate of global warming. The accurate estimation of forest carbon storage is related to the ability to reduce atmospheric CO₂ concentration. Against this background, the necessity of the study of forest carbon storage is revealed and its accurate estimation is an important basis for regional and even global carbon change studies.

The methods for estimating forest carbon storage includes both traditional direct measurement methods and indirect estimation methods based on remote sensing [3,4]. The traditional method is based on sample site survey and data statistics, which require a large investment of human and financial resources [5], and due to cost and time constraints, there are significant limitations regarding the large range of study objectives and long study periods. In contrast, remote sensing-based methods have the characteristics of being fast, low-cost, large-scale, and less destructive [1,6]. With the rapid development of quantitative remote sensing technology, the results of carbon storage estimation based on remote sensing have become more and more stable and reliable. Therefore, the use of remote sensing methods for the quantitative estimation of forest carbon storage on a larger scale has become a major current trend [7,8].

At present, there is a significant amount of research in remote sensing estimation, and while scholars in various fields have constructed different remote sensing models [9], most of these are static models. The static model is a model developed directly from the attribute values of the sample sites and the remotely sensed feature values of the corresponding points, which are uncalculated state data. If these values are further processed and calculated to produce different types of change values (periodic change, annual average change, rate of change, etc.), then the model developed using these change data is the dynamic model. Change data can describe the different processes of change in forest ecosystems at a given time and have the ability to monitor forest dynamics [10]. Gómez et al. [11] and Zhang et al. [12] compared static models with dynamic models and showed that dynamic models have higher accuracy and better predictive power than static models, with the highest R² for dynamic models in their study being 0.90 and 0.94, respectively. In the study of forest dynamics, the dynamic model is able to directly estimate the corresponding level of forest change and directly respond to forest dynamics, which improves the efficiency of the long-term forest dynamics study. The calculation of change data needs to be supported by multi-period data. The current National Forest Inventory (NFI) projects in many countries provide the possibility for long-term forest dynamics modelling studies. The ongoing NFI projects in many countries and open-access Landsat time series data offer the potential for long-term forest dynamics modelling studies.

Trees accumulate carbon mainly through photosynthesis during growth, and topographical factors influence the photosynthetic effect by affecting light and water conditions in the environment, which eventually affects the carbon storage of trees [13]. Therefore, topographic factors are an important concern in forest carbon storage studies. The digital elevation model (DEM) provides valuable topographic information including elevation, slope, and aspect [14]. These topographic factors are often used as independent factors in various studies because they do not require complex calculations and processing, are relatively easy to obtain, and reflect the effects of topography on different aspects of vegetation to some extent [15]. However, the influence of topography on forest carbon storage is the result of the combined effect of various topographic factors, and it is difficult to reflect this combined effect by a single elevation or slope. The terrain niche index (TNI) combines elevation and slope information, which is able to reflect the spatial variability of regional elevation and slope [16,17]. Therefore, in past studies, it has been frequently applied to land use change, landscape patterns and ecological effects [18,19]. The topography of forest ecosystems is complex and diverse, so the comprehensive influence of topographic factors should also be considered when modelling carbon storage. The change data can reflect the degree of change in forest carbon storage, and TNI can reflect the comprehensive influence of topographic factors on forest carbon storage, so the combination of the two is important for the accurate estimation of change in forest carbon storage in areas with complex topography.

Shangri-La, Yunnan Province, is rich in forest resources, and Pinus densata is one of the major tree species, which is widely distributed and has an essential impact on the carbon balance of the region. This region is located at the edge of the Tibetan Plateau, with an overall high altitude, and the terrain comprises high mountains and deep valleys with large topographic undulations, making it more difficult to estimate the forest carbon storage in this region accurately. Therefore, in this study, we used the remote sensing dynamic model combined with TNI to estimate the carbon storage of Pinus densata in Shangri-La. The main objectives are the following: (1) developing the carbon storage dynamics model based on two different kinds of change data; (2) exploring the impact of TNI on carbon storage models; and (3) estimating the historical carbon storage changes of Pinus densata in Shangri-La and analyzing the topographic effects of carbon storage changes.

2. Materials and Methods

2.1. Study Area and Study Process

The study area, Shangri-La, is located in Yunnan Province in southwestern China (Figure 1), at the border of the Yunnan Province and the Sichuan Province, and the geographical coordinates range from 99°20′ to 100°19′ eastern longitude and 26°52′ to 28°52′ northern latitude, the total area is 1.16 × 10⁶ ha. The area is located in the southeastern part of the Tibetan Plateau (the average elevation is 3459 m), with the Hengduan Mountains running north and south through the whole territory, the regional topography is undulating [20], and the whole terrain is high in the northwest and low in the southeast. Shangri-La has a remarkable monsoon climate, with rainfall concentrated in the months of June to October each year [21], the average annual precipitation is 268~945 mm, and the annual sunshine is 1742.9~2186.6 h. This region is very rich in forest resources, with 76% of forest cover [22]. Pinus densata, Pinus yunnanensis, and Picea asperata are the main tree species in the region [23]. The target species for this study is the Pinus densata, which includes both artificial and natural forests distributed within the study area.

The main content of the study is shown in Figure 2. The research process began with the pre-processing of all the data; then, two types of change data (annual average change and 5-year interval change) were calculated, three methods were used to construct the model, and TNI was added to the model. Finally, in order to compare the effect of the dynamic model with different change data and to compare the effect of the model before and after the addition of TNI, different indicators were used to evaluate the accuracy.

2.2. Ground Survey Data and Carbon Storage Calculation

This study focused on Pinus densata, and the ground survey data were obtained from the National Forestry Inventory, which contains a total of 136 Pinus densata sample plots (comprising pure forests of Pinus densata or forests with Pinus densata as the main species). The years of the collection were 1987, 1992, 1997, 2002, 2007, 2012 and 2017, where the number of sample plots was 19, 22, 23, 16, 16, 17 and 23 for each year, respectively (Figure 3). The dataset records information on tree height, diameter at breast height (DBH), number of trees, coordinate location, and major tree species. The distance between the sample plots consisted of regular distributions of 6 km × 8 km, and each sample plot was a rectangle of 28.28 m × 28.28 m (0.08 ha).

The aboveground biomass (AGB) of the sample plots was calculated based on the allometric growth equation of Pinus densata [24]. The average AGB was first calculated using the average tree height and average DBH of the sample plots. Then, the total AGB of the sample plots was calculated based on the average AGB and the number of Pinus densata. The allometric growth equation is as follows:

$${ \mathrm{AGB}}_1=0.073\times {\mathrm{DBH}}^{1.739}\times {\mathrm{H}}^{0.880}$$

(1)

where AGB₁ is single wood aboveground biomass (t), DBH is the diameter at breast height (cm), and H is tree height (m). We filtered the sample plot data. In this process, five sample plots with an AGB (AGB < 1 t·ha⁻¹) that was too small were removed, followed by six outliers that were screened out and removed using the Pauta criterion [25]. According to the Pauta criterion, if a value is outside the range of three times the standard deviation of the mean (outside the range of $\overline{\mathrm{x}}\pm 3\mathsf{\sigma }$, where $\overline{\mathrm{x}}$ is the mean and $\mathsf{\sigma }$ is the standard deviation), it is considered an outlier. This is a common method for outlier screening [26]. In the end, 125 sample plots of Pinus densata remained.

The carbon storage of the sample sites was calculated using the AGB multiplied by the carbon content coefficient. According to the guidelines for measuring carbon storage in forest ecosystems issued by the State Forestry and Grassland Administration of China [27], the average carbon content in the dry matter of Pinus densata is 0.501. The equation is as follows:

$$\mathrm{CS}=0.501\times {\mathrm{AGB}}_2$$

(2)

where CS is the carbon storage (t-C·ha⁻¹), and AGB₂ is the aboveground biomass of the sample plots (t·ha⁻¹). The basic sample data is shown in Figure 4. The box plots (Figure 4) depict information on the diameter at breast height, tree height, and carbon storage for each year of the sample plots. The overall range of DBH is 6.0–91.7 cm, with a maximum annual mean of 32.0 cm and a maximum annual standard deviation of 27.8 cm; the overall range of tree height is 1.5–24.8 m, with a maximum annual mean of 11.4 m and a maximum annual standard deviation of 5.5 m; and the overall range of carbon storage is 0.7–84.2 t-C·ha⁻¹, with a maximum annual mean of 38.6 t-C·ha⁻¹ and a maximum annual standard deviation of 41.9 t-C·ha⁻¹.

2.3. Remote Sensing Images and Obtained Features

2.3.1. Remote Sensing Data

The remote sensing images used in this study were obtained from the USGS website (http://glovis.usgs.gov/ (accessed on 28 October 2021)). Landsat 5 TM images from 1987, 1992, 1997, 2002, 2007, and 2012 and Landsat 8 OLI images from 2017 of the Shangri-La region were downloaded from the website (

Table 1. Landsat time-series of imagery used in the study.

Year	Landsat/Sensor	Path/Row	Image Acquisition Date
1987	5/TM	132/040	30 December 1987
	5/TM	132/041	30 December 1987
	5/TM	131/041	23 December 1987
1992	5/TM	132/041	7 November 1991
	5/TM	132/040	7 November 1991
	5/TM	131/041	16 November 1991
1997	5/TM	132/041	7 November 1997
	5/TM	132/040	6 October 1997
	5/TM	131/041	16 November 1997
2002	5/TM	132/041	5 January 2002
	5/TM	132/040	5 January 2002
	5/TM	131/041	29 October 2002
2007	5/TM	132/041	15 October 2006
	5/TM	132/040	3 January 2007
	5/TM	131/041	1 March 2007
2012	5/TM	132/041	13 October 2011
	5/TM	132/040	14 January 2011
	5/TM	131/041	7 January 2011
2017	8/OLI	132/041	16 December 2017
	8/OLI	132/040	16 December 2017
	8/OLI	131/041	25 December 2017

), comprising 21 views in total and a spatial resolution of 30 m. When an image contains too much cloud or is of poor quality, we chose an image from a neighboring timepoint to replace it. Most of the selected images contained less than 5% cloud, and the cloud did not cover the study area when the cloud content was greater than 5%.

In order to improve the image quality, all the images were preprocessed: firstly, radiometric calibration was performed, which converts the original DN (digital number) value of the image to a consistent radiometric brightness in order to eliminate the influence of the sensor [28]; secondly, the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) method was used to eliminate the influence of the atmosphere [29]; thirdly, geometric correction was performed on each image with reference to the existing standard SPOT-5 images to eliminate geometric errors, where the error in geometric correction was less than one pixel; and finally, in order to eliminate the change of image grey value caused by the topographic undulation, topographic correction was performed using the slope matching model [30]. After topographic correction, the brightness of the illuminated side was suppressed and the brightness of the shadow side was enhanced, which better restored the reflectance of the hidden features in the shadow of the original image [31].

2.3.2. Remote Sensing Features

The remote sensing pixel values were evaluated by considering bands with the same name and close wavelengths in Landsat 8 and Landsat 5. In this study, two types of remote sensing features, spectral features and texture features, are calculated. Spectral features are the most basic and direct features of remote sensing images, but the phenomena of “the same objects with different spectra” and “different objects with the same spectra” commonly exist in remote sensing images, so it may be difficult to obtain good research results when only relying on spectral features [32,33]. The texture features can describe the grey scale information of an image, express the spatial distribution of the grey scale of pixels in the image, and reflect the structural information of the image. Compared with spectral features, texture features are less affected by the environment and more stable in information representation [34]. Therefore, the combination of spectral and textural features can effectively reflect the characteristics of the ground objects and their changes [35,36]. A total of 35 spectral features and 540 texture features were obtained in this study (the calculations of various types of textures from all single bands of each image in odd windows from 1 to 19) (

Table 2. Extracted remote sensing feature information.

Categories		Information on Remote Sensing Characteristics
Spectral features	Original bands	B1; B2; B3; B4; B5; B7
	Vegetation indices [37]	$\mathrm{NDVI}=\left(\mathrm{B}4-\mathrm{B}3\right)/\left(\mathrm{B}4+\mathrm{B}3\right)$; $\mathrm{ND}32=\left(\mathrm{B}3-\mathrm{B}2\right)/\left(\mathrm{B}3+\mathrm{B}2\right)$; $\mathrm{ND}54=\left(\mathrm{B}5-\mathrm{B}4\right)/\left(\mathrm{B}5+\mathrm{B}4\right)$; $\mathrm{ND}53=\left(\mathrm{B}5-\mathrm{B}3\right)/\left(\mathrm{B}5+\mathrm{B}3\right)$; $\mathrm{ND}57=\left(\mathrm{B}5-\mathrm{B}7\right)/\left(\mathrm{B}5+\mathrm{B}7\right)$; $\mathrm{ND}452=\left(\mathrm{B}4+\mathrm{B}5-\mathrm{B}2\right)/\left(\mathrm{B}4+\mathrm{B}5+\mathrm{B}2\right)$; $\mathrm{DVI}=\mathrm{B}4-\mathrm{B}3$; $\mathrm{RVI}=\mathrm{B}4/\mathrm{B}3$; $\mathrm{ARVI}=\left(\mathrm{B}4-\left(2\mathrm{B}3-\mathrm{B}1\right)\right)/\left(\mathrm{B}4+\left(2\mathrm{B}3-\mathrm{B}1\right)\right)$ $\mathrm{EVI}=2.5\times \frac{\mathrm{B}4-\mathrm{B}3}{\mathrm{B}4+6\mathrm{B}3-7.5\mathrm{B}1+1}$;
	Band combination [12]	$\mathrm{B}4/\mathrm{B}2$; $\mathrm{B}5/\mathrm{B}3$; $\mathrm{B}5/\mathrm{B}4$; $\mathrm{B}5/\mathrm{B}7$; $\mathrm{B}7/\mathrm{B}3$; $\mathrm{B}4\times \left(\mathrm{B}3/\mathrm{B}7\right)$
	Image enhancement [38,39]	$\mathrm{Albedo}=\mathrm{B}1+\mathrm{B}2+\mathrm{B}3+\mathrm{B}4+\mathrm{B}5+\mathrm{B}7$; $\mathrm{VIS}123=\mathrm{B}1+\mathrm{B}2+\mathrm{B}3$; $\mathrm{MID}=\mathrm{B}5+\mathrm{B}7$; $\mathrm{FVC}=\frac{\mathrm{EVI}-{\mathrm{EVI}}_{\min }}{{\mathrm{EVI}}_{\max }-{\mathrm{EVI}}_{\min }}$; Principal component analysis (PCA1, PCA2, PCA3, PCA4, PCA5, PCA7); Tasseled cap transform (TCT1, TCT2, TCT3)
Texture features	Grey level co-occurrence matrix [34]	Homogeneity (HO); Dissimilarity (DI); Mean (ME); Angular second moment (SM); Entropy (EN); Correlation (CC); Variance (VA); Contrast (CO)
	Filtering of probabilistic statistics [40]	Skewness (SK)

Where B1 is the blue band, B2 is the green band, B3 is the red band, B4 is the near-infrared band, B5 is the shortwave infrared-1 band, and B7 is the shortwave infrared-2 band.

).

2.4. Terrain Niche Index

To explore the influence of elevation and slope on the carbon storage dynamics model, the terrain niche index (TNI), a composite factor combining elevation and slope, was established in this study [41]. Elevation and slope were extracted from the digital elevation model (DEM). The DEM used in this study was an ASTER GDEM data product with a spatial resolution of 30 m. The DEM was downloaded from the Geospatial Data Cloud website (http://www.gscloud.cn/ (accessed on 18 August 2022)). The equation for TNI is as follows:

$$\mathrm{TNI}=\log \left[\left(\frac{\mathrm{E}}{\overline{\mathrm{E}}}+1\right)\times \left(\frac{\mathrm{S}}{\overline{\mathrm{S}}}+1\right)\right]$$

(3)

where $\mathrm{E}$ and $\mathrm{S}$ are the elevation and slope of any point, respectively, and $\overline{\mathrm{E}}$ and $\overline{\mathrm{S}}$ are the average elevation and average slope of the study area, respectively. A larger TNI value indicates a higher elevation and slope; a smaller TNI value indicates a lower elevation and gentler slope; and a medium TNI value indicates a high elevation gentle slope, or a low elevation steep slope, or a medium elevation and slope [42].

2.5. Distribution Index

The distribution index (DI) is often used to show the distribution of various land use types. In this study, the areas of the increasing and decreasing carbon storage of Pinus densata were considered as two different land classes. Then, the DI was used to explore the distribution of carbon storage changes in different TNI gradients. The DI is a standardized, dimensionless metric that eliminates the effects of area differences and allows the comparison of distribution characteristics between carbon storage changes of different area proportions [18,43]. The equation is as follows:

$$\mathrm{P}=\left({\mathrm{S}}_{\mathrm{ie}}/{\mathrm{S}}_{\mathrm{i}}\right)/\left({\mathrm{S}}_{\mathrm{e}}/\mathrm{S}\right)$$

(4)

where $\mathrm{P}$ is the DI, $\mathrm{e}$ is the gradient of TNI, ${\mathrm{S}}_{\mathrm{ie}}$ is the area of $\mathrm{i}$ changes on the $\mathrm{e}$ TNI gradient, ${\mathrm{S}}_{\mathrm{i}}$ is the area of the $\mathrm{i}$ change, ${\mathrm{S}}_{\mathrm{e}}$ is the area of the TNI for gradient $\mathrm{e}$, and $\mathrm{S}$ is the total area of the study area. A larger $\mathrm{P}$ indicates a higher frequency of this type of change, and $\mathrm{P}>1$ indicates that the type of change under this TNI gradient belongs to the dominant distribution. Therefore, the analysis of the distribution of each carbon storage change type on different TNI levels is able to reveal the influence of different topographic conditions on carbon storage change.

2.6. Modelling Process

2.6.1. Calculation of Change Data

In this study, two kinds of change data were used to develop the forest carbon storage model, one is the 5-year interval change, and the other is the annual average change. The change data were calculated based on a common sample plot for both years: we calculate the change value of a sample plot when data were recorded in both adjacent years, otherwise the change value for the sample plot during that specific time period is not calculated. The equations are as follows:

$$∆\mathrm{I}={\mathrm{I}}_{\mathrm{n}}-{\mathrm{I}}_{\mathrm{m}}$$

(5)

$$∆\mathrm{R}=\frac{∆\mathrm{I}}{\mathrm{n}-\mathrm{m}}$$

(6)

where $∆\mathrm{I}$ is the interval change value, $\mathrm{n}$ and $\mathrm{m}$ are two different years in which the change is calculated ($\mathrm{n}>\mathrm{m}$), ${\mathrm{I}}_{\mathrm{n}}$ and ${\mathrm{I}}_{\mathrm{m}}$ are the data values for years $\mathrm{n}$ and $\mathrm{m}$, respectively, and $∆\mathrm{R}$ is the average annual change value. The above two equations were used to calculate the change data of carbon storage and the corresponding remote sensing features in the sample sites. The calculation of change values is supported by the NFI data records every 5 years from 1987–2017. The shortest time interval for which change data can be calculated is 5 years and the longest time interval is 30 years. In generally, the shorter the time interval, the stronger the continuity of the data and the more accurately the dynamic process of forest change can be expressed [10,44]. Zhang et al. [12] tested the change data for 5, 10, and 15 year intervals and found that the 5 year interval resulted in the smallest error. Lunetta et al. [45] compared the effect of 3, 5, and 7 year intervals on monitoring change and showed that the 3 year interval had the highest accuracy. Therefore, only the change data for the shortest interval (5 years) were calculated in this study: the 5-year interval change and the annual average change during the 5-year period.

2.6.2. Remote Sensing Features Selection

The efficiency of the model can be affected by the number of features, because too many features reduce the speed of model fitting [46]. In addition, the original features contain some redundant information, which has a negative impact on the model accuracy [47]. Therefore, the feature screening was performed. In this study, 575 extracted remote sensing features were screened using the feature importance assessment method of the random forest [48]. In order to show the combined effect of spectral features and texture features in modelling, the spectral feature and texture feature were screened independently in this study. These two features were then combined and used for modelling.

Since most of the texture values of the extracted 1 × 1 windows were 0 or 1, they were not usable and thus removed, and the remaining 486 texture features were used for the subsequent screening. The 5-year interval change data and the annual average change data were calculated using Equations (5) and (6), respectively. Then, the feature importance assessment method of the random forest was used to filter the feature variables in the two different change types, and the top five features with contributions greater than 5% were selected from the spectral features and the texture features, respectively. Considering the problem of collinearity between features, the variance inflation factor (VIF) was used to check all the features selected in this study. When the VIF value is less than 3, it means that the collinearity between features is weak [49]. According to this principle, two features with VIF values greater than 3 in the 5-year interval change were removed, and the 5-year interval change and annual average change resulted in 8 and 10 features for the modelling study, respectively (

Table 3. Results of remote sensing feature screening.

Change Type	Remote Sensing Features
5-year interval change	$\mathrm{B}4\times \left(\mathrm{B}3/\mathrm{B}7\right)$, $\mathrm{B}5/\mathrm{B}4$, PCA2, FVC, PCA4, R9B1EN, R17B7CO, R5B4SM
Annual average change	R3B3VA, R15B5VA, R17B1VA, R9B4HO, R15B5ME, DVI, PCA2, EVI, $\mathrm{B}4\times \left(\mathrm{B}3/\mathrm{B}7\right)$, $\mathrm{B}7/\mathrm{B}3$

The expression of the texture features is “RXBYTT”, RX is the size of the texture window, BY is a certain single band of the image, and TT is the abbreviation of a certain texture feature.

).

2.6.3. Modelling Methods

In this study, carbon storage dynamics models were developed using change data, and the modelling methods included partial least squares regression (PLSR), gradient boosting regression tree (GBRT), and random forest (RF).

PLSR is a method used to study the correlation and quantitative relationship between the response variable $Y$ and a set of explanatory variables, $X=x_1,x_2,\cdots ,x_n$, and the set of independent variables $X$ can also be used to predict $Y$ [50]. The modelling principle of PLSR is the development of multivariate linear regression and principal component analysis. Compared with these two methods, PLSR has better stability and can solve the problem of collinearity between multiple explanatory variables $X$ [51,52,53]. The implementation of the PLSR model in this study was based on the Minitab 20 software [54].

GBRT is an algorithm for iterative regression trees proposed by Friedman [55,56]. All the regression trees in this model are interconnected, and at its core, each tree is fitted based on the residuals and conclusions of the previous tree [57]. GBRT uses a forward-distributed algorithm to minimize the loss function by selecting the appropriate decision tree function from the current model and the fitted function. Based on the characteristics described above, the GBRT model has low computational complexity, is able to reduce errors, and has the ability to handle unevenly distributed data [58]. The implementation of the GBRT model in this study is based on the “Gradient Boosting Regressor” algorithm provided in the “Scikit-learn” package for the Python language.

RF is an integrated learning model, first proposed by Breiman, which consists of many aggregated randomly generated trees [59,60]. In the classification problem, RF outputs the type with the most votes, and in the regression problem, RF outputs the mean of all decision trees [61]. RF can reduce variance and effectively reduce overfitting by assembling different trees, and it has excellent classification and regression performance; thus, it is currently used in many fields of research [58,61,62]. The implementation of the RF model in this study is based on the “Random Forest Regressor” algorithm provided in the “Scikit-learn” package for the Python language.

In this study, the above three methods were used to develop 5-year interval change and annual average change models, respectively. In order to explore the influence of TNI on the dynamic model of forest carbon storage, the TNI was added to each model in this study, then the accuracy of different dynamic models before and after adding the TNI was compared, and the optimal model was selected to estimate the carbon storage of Pinus densata in Shangri-La.

2.7. Accuracy Evaluation

Ultimately, this study obtained 92 groups of change values of each type of change; 70% (64 groups) of the data were randomly selected for model fitting, and the remaining 30% (28 groups) were used for validation, and cross-validation was performed during model fitting. The indicators used to evaluate the accuracy of the model were the coefficient of determination (R²), the root mean square error (RMSE), and the mean absolute error (MAE). In order to ensure that the model results were as objective as possible, each model was fitted 20 times in this study to allow take the mean values of evaluation indicators to be used for comparison. The equations are as follows:

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

(7)

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

(8)

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{y}}_{\mathrm{i}}-{\hat{\mathrm{y}}}_{\mathrm{i}}\right|$$

(9)

where ${\mathrm{y}}_{\mathrm{i}}$ is the observed value of the sample site, ${\hat{\mathrm{y}}}_{\mathrm{i}}$ is the predicted value of the model, $\overline{\mathrm{y}}$ is the observed mean value of the sample site, and n is the number of samples.

3. Results

3.1. Analysis of Modelling Results

The PLSR, GBRT, and RF methods were used to develop the 5-year interval change and annual average change models of carbon storage, respectively, TNI was added to each model subsequently, and the results of accuracy evaluation indexes (20 times fitted mean values) of each type of model are shown in

Table 4. Comparison of modelling results.

Change Type	Model	Fitting		Validation
		R2	RMSE/t-C·ha−1	MAE/t-C·ha−1
5-year interval change	PLSR	0.18	9.61	9.09
	GBRT	0.81	4.58	2.51
	RF	0.85	4.09	1.83
	PLSRTNI	0.18	9.50	8.36
	GBRTTNI	0.83	4.24	2.43
	RFTNI	0.87	3.82	1.78
Annual average change	PLSR	0.23	1.83	1.64
	GBRT	0.80	0.94	0.45
	RF	0.83	0.87	0.36
	PLSRTNI	0.25	1.82	1.67
	GBRTTNI	0.82	0.89	0.41
	RFTNI	0.84	0.86	0.35

PLSR _TNI, GBRT_TNI, and RF_TNI in the table indicate the PLSR model, GBRT model, and RF model after adding TNI, respectively.

. As can be seen from Table 4, among the models developed with the same type of change data, the PLSR model performed the worst and the RF model performed the best. The change of model indicators (percentage) after adding TNI is shown in Figure 5, from which it can be seen that R² increases and RMSE decreases for all models, and the MAE values decreased for all models except for the PLSR model with annual average change. The maximum increase in R² is 10.48%, the maximum decrease in RMSE is 7.32%, and the maximum decrease in MAE is 8.89%. Considering the changes of all accuracy indicators together, the addition of TNI improves the model accuracy. The mean values of carbon storage calculated for the 5-year interval change and the mean annual change of the sample plot data were 5.14 t-C·ha⁻¹ and 1.03 t-C·ha⁻¹, respectively. After comparison, it can be seen that the RMSE and MAE values of both the RF and GBRT models are smaller than the mean values, except for the PLSR model, indicating that the errors of the RF and GBRT models are smaller and the model results are more reasonable.

The 5-year interval change represents the change in carbon storage over five years, and the annual average change represents the change in carbon storage in one year. In order to compare the effects of dynamic models with 5-year interval changes and annual average changes, the RMSE and MAE values of the annual average change model were adjusted to five times their original values in this study, so that the error indicators of both models are on the same time scale. The adjusted results are shown in

Table 5. Error indicator for the adjusted annual average change.

Change Type	Model	Fitting	Validation
		RMSE/t-C·ha−1	MAE /t-C·ha−1
Annual average change	PLSR	9.16	8.2
	GBRT	4.69	2.25
	RF	4.37	1.8
	PLSRTNI	9.075	8.35
	GBRTTNI	4.465	2.05
	RFTNI	4.295	1.75

. Comparing the accuracy indicators of the two models in Table 4 and Table 5, it can be seen that the R² of all 5-year interval change models is higher than that of the annual average change model, the RMSE values are lower than those of the 5-year interval change model except for the PLSR model, and the difference in MAE values is not significant. Therefore, the 5-year interval change model works better overall.

According to the confidence interval plot (Figure 6) of the 20 fittings for each model, the confidence range of the RF model among the types is smaller than that of the GBRT and PLSR models, indicating that the RF model in this study has the least uncertainty and the results obtained are more stable and reliable. In summary, the accuracy of 5-year interval change RF model with the addition of TNI is the highest.

3.2. Mapping Carbon Storage

According to the model evaluation results, the 5-year interval change RF model with the highest accuracy was used to estimate the carbon storage change. The scatter plot of this estimation model is shown in Figure 7. The direct estimation result of the model is the 5-year interval change of carbon storage of Pinus densata in Shangri-La. The values of the carbon storage change for the six change periods are ${\mathsf{\Delta }\mathrm{CS}}_1={\mathrm{CS}}_{1992}-{\mathrm{CS}}_{1987}$, ${\mathsf{\Delta }\mathrm{CS}}_2={\mathrm{CS}}_{1997}-{\mathrm{CS}}_{1992}$, ${\mathsf{\Delta }\mathrm{CS}}_3={\mathrm{CS}}_{2002}-{\mathrm{CS}}_{1997}$, ${\mathsf{\Delta }\mathrm{CS}}_4={\mathrm{CS}}_{2007}-{\mathrm{CS}}_{2002}$, ${\mathsf{\Delta }\mathrm{CS}}_5={\mathrm{CS}}_{2012}-{\mathrm{CS}}_{2007}$ and ${\mathsf{\Delta }\mathrm{CS}}_6={\mathrm{CS}}_{2017}-{\mathrm{CS}}_{2012}$.

The estimated carbon storage changes in the six periods are as follows: ΔCS₁ = −53.327 × 10⁴ t, ΔCS₂ = 13.887 × 10⁴ t, ΔCS₃ = −94.041 × 10⁴ t, ΔCS₄ = 14.602 × 10⁴ t, ΔCS₅ = 20.869 × 10⁴ t, and ΔCS₆ = 5.602 × 10⁴ t. A positive value indicates an increase in total carbon storage, and a negative value indicates a decrease in total carbon storage. The estimation results indicate that the carbon storage of Pinus densata in Shangri-La decreased in the periods 1987–1992 and 1997–2002 and increased in the remaining four periods.

In order to obtain carbon storage maps for each year, the published AGB estimates for each year from our group were used [63], which contain the aboveground biomass estimates of Pinus densata in Shangri-La for 1987, 1992, 1997, 2002, 2007, 2012, and 2017. These data were combined with the results of carbon storage change values estimated in this paper to calculate the total carbon storage values for each year. Firstly, the AGB of each year in the published study was converted to carbon storage by multiplying the carbon content coefficient (0.501) and the converted values were expressed as T₁₉₈₇, T₁₉₉₂, T₁₉₉₇, T₂₀₀₂, T₂₀₀₇, T₂₀₁₂, and T₂₀₁₇; then the carbon storage change values ΔCS₁, ΔCS₂, ΔCS₃, ΔCS₄, ΔCS₅ and ΔCS₆ were added to the conversion values of the smaller years in each change period in order to obtain the carbon storage values for the six years 1992, 1997, 2002, 2007, 2012 and 2017, respectively (e.g., ${\mathrm{CS}}_{1992}={\mathsf{\Delta }\mathrm{CS}}_1+{ \mathrm{T}}_{1987}$), and the carbon storage value for 1987 was obtained by subtracting ΔCS₁ from T₁₉₉₂ (${\mathrm{CS}}_{1992}={\mathrm{T}}_{1987}-{\mathsf{\Delta }\mathrm{CS}}_1$). The statistical results for carbon storage are shown in

Table 6. Statistics of Pinus densata carbon storage.

Year	Area of Pinus densata (ha)	Total CS (Million Tons)	Average CS (t-C·ha−1)
1987	171560.28	5.30	30.91
1992	171560.28	4.24	24.69
1997	170589.86	4.77	27.97
2002	170589.86	3.12	18.30
2007	174179.37	3.99	22.93
2012	174213.12	4.03	23.12
2017	184815.84	3.80	20.53

, and Figure 8 and Figure 9 (the spatial resolution is 30 m) were obtained by mapping the values of the carbon storage change for each period and each year according to the Pinus densata distribution range. The distribution range and area of Pinus densata are derived from Forest Manager Inventory (FMI) data.

3.3. Spatial Distribution Characteristics for the Changes in Carbon Storage on Different TNI Gradients

In this study, the estimated changes in carbon storage were classified as increases (positive change values) and decreases (negative change values), and the TNI (0.30 to 0.93) was classified into ten levels from low to high using the natural breaks (Jenks) method: low (1–3), medium–low (4–5), medium–high (6–7), high (8–10). The DI was then calculated according to Equation (4) and used to explore the distribution of carbon storage changes along different TNI gradients. The distribution results are shown in Figure 10.

As can be seen from Figure 10, the increase and decrease in the carbon storage of Pinus densata in Shangri-La at each period have a certain regularity on the TNI gradient. The DI curve for ΔCS₁(1987–1992) decreases and then rises when carbon storage increases (Figure 10a) with a dominant distribution in the low and high TNI gradients; the DI curve for ΔCS₂ (1992–1997), ΔCS₄ (2002–2007), ΔCS₅ (2007–2012), and ΔCS₆ (2012–2017) rises and then reduces with the dominant distribution areas of ΔCS₂ being in the low and medium TNI gradients, and the dominant distribution areas of ΔCS₄, ΔCS₅, and ΔCS₆ being concentrated in the medium–high TNI gradients; and ΔCS₃ (1997–2002) reduces with the TNI gradient, and its dominant distribution area is in the low and medium–low gradients. The corresponding curve change period displays a decreasing trend when the carbon storage decreases (Figure 10b). Overall, the dominant distribution areas for each period when carbon storage increases are mainly found in the lower or middle TNI gradients, while the dominant distribution areas for each period when carbon storage decreases are found in all TNI gradients.

4. Discussion

In long time series data, the uncertainty of the historical data, due to the events that took place at the time being taken into account, cause some of the data to display large deviations in values, which are known as anomalies or outliers [64]. The stability of the model can be affected by outliers, and at the same time, the prediction accuracy of the model may be reduced, so the detection and handling of outliers is very important in the modelling process [65]. The sample plot data were screened twice in this study, the first time removing plots with an AGB of less than 1 t·ha⁻¹ because they had an average diameter at breast height of less than 5 cm and an average tree height of less than 1.5 m. They are young forests, too low in canopy density, and poorly characterized by the appearance of the forest floor (similar to bare ground) on 30 m resolution imagery. Values exceeding the range of $\overline{x}\pm 3\sigma$ were selected for the second time using the Pauta criterion. By removing values that are too small and too large from the data in these two steps, respectively, the overall uncertainty of the data is reduced and the reliability of the data is improved.

The remote sensing-based estimation model is one of the main methods used in current carbon storage studies [7]. Based on the accumulation of ground survey data, the results of remote sensing estimation models are considered to have certain advantages and reliability in forest carbon storage studies. Time series images also provide an important basis for describing change [66]. The NFI and Landsat time series data over a 30-year period were combined to develop a dynamic model to study forest carbon storage in this paper. Although dynamic models based on change data are currently less commonly used in forestry, their advantages in terms of accuracy have been demonstrated [11,12]. In this study, although the effects of dynamic models and static models are not directly compared through experiments, we can refer to some of the existing studies related to forest carbon storage modelling in recent years [8,67,68]. The R² of the remote sensing static models developed in these studies ranges from 0.64 to 0.73 as the highest value, while the R² of the optimal model in this study is 0.87. Thus, the advantages of the dynamic model in terms of accuracy can be seen.

The DI is an area evaluation indicator. When the distribution curve is flatter, it indicates that the distribution of such changes deviates less from the standard distribution, and its adaptability to the terrain is wider [16]. In this study, the areas of increasing or decreasing Pinus densata carbon storages were calculated separately, and the main topography of these two areas was evaluated according to the DI. The results of this study demonstrate that the increase in carbon storage in the periods ΔCS₄, ΔCS₅, and ΔCS₆ is much greater than the decrease, leading to an increase in total carbon storage from 2002 to 2017, which is related to the policy of “returning farmland to the forest” that has been implemented in Shangri-La since 2000. In the area of increasing carbon storage, the distribution curves of these three periods are relatively flat, showing a general adaptation to different levels of TNI and suggesting that the implementation of this government policy has led to the widespread planting of Pinus densata in Shangri-La, thus resulting in a more balanced topographic distribution of increasing carbon storage across the region.

Different TNI levels reflect different elevation and slope conditions, and the dominant distribution reflects the main topographic range of carbon storage changes. Although the dominant distribution of carbon storage changes over the TNI gradient can be clearly seen throughout a particular period of change, the six periods of change in carbon storage are more complex, and the distribution trends across the periods of change are not entirely consistent along the TNI gradient. What can be seen from the dominant distribution is that the regions where carbon storage increased are mainly located at low altitudes or on gently sloping terrain or both, as such terrain is conducive to artificial tending. In contrast, in regions with decreased carbon storage, the dominant distribution occurs across the TNI gradient, and the overall dominance of the terrain is not obvious, which is due to the randomicity of deforestation or forest destruction. In the past, TNI has not been applied to separate forest land class studies, but this study shows that its combination with the DI can be used to evaluate the increase and decrease effects of a single land type.

Some studies have suggested that human activities are the main driver of forest carbon [69], but as Shangri-La is a highland region with a low population density, the impact of human life on the surrounding forests is also relatively low, and it is environmental factors, including topography, that mainly impact forest carbon storage. Therefore, it is essential to consider topographic factors when studying forest carbon storage in this region. The results of this paper show that adding TNI to the dynamic model can reflect the topographic effect of the model and enable the improvement of the dynamic model accuracy, and by analyzing the distribution of forest carbon storage changes along the TNI gradient, it is possible to reveal the changes in forest carbon storage under different topographic conditions.

At present, the main approach to the study of historical forest dynamics is still based on the estimation of forest biomass or carbon storage over several single years, and the values of changes in different periods are calculated in order to derive forest dynamics results [70]. The level of periodic change in the forest can be estimated directly from models developed from interval change data, and due to the advantages of dynamic models in terms of accuracy, the resulting change values are estimated more precisely. Different dynamic models can be obtained for different types of change data. In this paper, we only compared two dynamic models for the 5-year interval change data and annual average change data, and more types of dynamic models are yet to be explored.

Although it was difficult to derive carbon storage results for a single year in this study without the support of carbon storage inventory data or graphs related to the study years, the results of this study demonstrate that the dynamic model is able to estimate the change in forest carbon storage and is well suited for the study of forest carbon storage change. Dynamic models are capable of quantitatively describing changes in specific properties, so they need not be limited to the field of forestry research but have the potential for application in other areas of natural or manufactured landscapes where regular variation exists. This study has so far only considered the influence of topographic factors, and there remains the possibility of other environmental factors, such as climate and soil, also influencing the model, which could be the next direction for future research.

5. Conclusions

In this study, the carbon storage dynamic model was developed based on NFI data and Landsat time series images over a 30-year period. PLSR, RF, and GBRT were used for modelling, and the accuracy of the two dynamic models was compared; then TNI, which represents topographic factors, was introduced into the model, and the distribution of carbon storage changes in Shangri-La on different TNI gradients was analyzed. The main conclusions are as follows: (1) the model effects of the non-parametric methods, RF and GBRT, are much better than those of the parametric method, PLSR, and the accuracy of the 5-year interval change model is better than that of the annual average change model; (2) TNI can improve the accuracy of the dynamic carbon storage estimation model, and the accuracy of the dynamic model with RF is the highest after adding TNI; (3) the dynamic model displays good performance regarding the estimation of carbon storage changes, and the results of the interval change model show that the total carbon storage of Pinus densata in Shangri-La decreased in 1987–1992 and 1997–2002, and increased in 1992–1997, 2002–2007, 2007–2012, and 2012–2017; and (4) the DI can be used to evaluate the main topography for the regions that display a change in Pinus densata carbon storage, the predominant topography in the regions of increasing carbon storage is that of low elevations or low slopes, or a combination of both conditions, while the topography in regions of decreasing carbon storage is more random. The results of this study can be used as a reference for forest carbon storage estimation using Landsat images in areas with complex topography.