Simulation of Spatial and Temporal Distribution of Forest Carbon Stocks in Long Time Series—Based on Remote Sensing and Deep Learning

Abstract

Due to the complexity and difficulty of forest resource ground surveys, remote-sensing-based methods to assess forest resources and effectively plan management measures are particularly important, as they provide effective means to explore changes in forest resources over long time periods. The objective of this study was to monitor the spatiotemporal trends of the wood carbon stocks of the standing forests in the southeastern Xiaoxinganling Mountains by using Landsat remote sensing data collected between 1989 and 2021. Various remote sensing indicators for predicting carbon stocks were constructed based on the Google Earth Engine (GEE) platform. We initially used a multiple linear regression model, a deep neural network model and a convolutional neural network model for exploring the spatiotemporal trends in carbon stocks. Finally, we chose the convolutional neural network model because it provided more robust predictions on the carbon stock on a pixel-by-pixel basis and hence mapping the spatial distribution of this variable. Savitzky–Golay filter smoothing was applied to the predicted annual average carbon stock to observe the overall trend, and a spatial autocorrelation analysis was conducted. Sen’s slope and the Mann–Kendall statistical test were used to monitor the spatial trends of the carbon stocks. It was found that 59.5% of the area showed an increasing trend, while 40.5% of the area showed a decreasing trend over the past 33 years, and the future trend of carbon stock development was plotted by combining the results with the Hurst exponent.

1. Introduction

Forests are the largest organic carbon pool in terrestrial ecosystems [1]. Although forests cover only one-third of the total land area, the carbon stock in forested areas accounts for 56% of the total terrestrial carbon pool [2]. In recent years, the sustainable management of forests has become an important tool to combat climate change, and changes in forest resources related to climate change and anthropogenic disturbances also affect global climate change [3]. Therefore, the estimation of forest carbon stocks is of critical importance for understanding the global carbon cycle and climate change. Forests offer great potential for the sequestration of atmospheric carbon, and thus, there has been an urgent need to accurately and efficiently quantify forest carbon stocks in recent years [4]. However, it is difficult to conduct field-based surveys of forest biomass and carbon stocks, and frequent field visits can further damage forest ecosystems. Therefore, it is challenging to estimate forest carbon sequestration potential or assess its spatial and temporal variability [5].

The emergence of remote sensing technology has provided additional possibilities for monitoring changes in forest resources, and its advantages, including rapid, real-time and large-scale monitoring, have made remote sensing a popular technique that is widely used in the fields of ecology and environment monitoring [6]. A series of remote sensing indices, such as NDVI [7] and EVI [8], used to describe and estimate the changes in forest resources, have provided a breakthrough in what was originally a relatively difficult problem. In addition, with the improvements in space technology, remote sensing indices such as satellite-derived land surface temperature (LST) [9,10] and wetness (WET) [11] have been proposed to effectively describe the changes in ecological indicators such as surface temperature and humidity. To find a suitable remote sensing factor to estimate forest carbon stock, we use Google Earth Engine (GEE) [12] to collect and process the required dataset. Additionally, GEE was used to perform the initial correction, remove clouds and provide a multitemporal synthesis of remote sensing images. Furthermore, GEE can convert the raw remote sensing images into our required remote sensing metrics. As a nonprofit open access platform, GEE [13] has been widely used in scientific research and education, especially in large-scale and long time series ecological evaluation, where it demonstrates diverse datasets and powerful computing capabilities [14].

The study site is located at the southeastern end of this ecosystem [15] and is a typical forest representative of the Xiaoxinganling Mountains that has strong carbon sequestration potential. Additionally, the reserve has both pristine forest, meadow and swamp ecosystems, as well as secondary, heavily disturbed and artificial forest ecosystems, which also have a richly diverse ecosystem and forest types [16]. Therefore, the reserve is already a mature natural experimental system, and relevant biomass studies have already been carried out in the region, but only for single forest types or biomass estimates for short time periods. The richness of forest types in the reserve and long-term quantitative surveys provide opportunities to study long time series forest carbon stock changes. Such comprehensive studies are essential for carbon stock assessment and effective forest management in the study area and throughout the Xiaoxinganling Mountains [17].

The emergence of a deep learning algorithm, an analysis technique, has further improved the usefulness of remote sensing data. Deep learning (DL) [18] represents a new research direction in the field of machine learning (ML) [19], which consists of neurons that gradually derive initial low-level feature representations into high-level feature representations through multilayer processing while steadily learning higher-level features [20]. Compared with statistical models and traditional machine learning techniques, DL is more capable of portraying and mining the rich intrinsic information contained within remote sensing data by transforming the feature representation of samples in the original space to a new feature space via layer-by-layer feature transformation [21]. Thus, the model uses remote sensing data to learn features.

The slope estimation is generally performed using the least squares method, and then the trend of the time series carbon stock is judged based on the slope [22]. However, due to the randomness and discontinuity of the data and the existence of outliers, the estimated slope does not usually detect that trend well. Sen’s slope estimator was proposed in 1968 for the nonparametric estimation algorithm of linear trends [23,24]. The advent of Sen’s slope estimator avoids the influence of outliers and distribution patterns in time series data on the analysis results. The Mann–Kendall (MK) statistical test is also used to determine the significance of trends in the time series [25]. MK is a nonparametric statistical test that does not require neither the normal distribution of the measurements nor a linear trend in the data, and it is not affected by missing values and outliers [26]. Sen’s slope estimator and the MK statistical test are widely used in hydrological and meteorological time series analysis, but their application in forest resource trend analysis is rare [27]. The Hurst exponent [28] is an analytical method for determining the long-term memory of time series. That work proposed the rescaled polar difference (R/S) analysis methods to establish the Hurst exponent (H) as an indicator to determine whether the time series data follow a random walk or a biased walk trend, and it is now widely used in ecological assessments [29].

In this context, we conducted a study to investigate the forest carbon stocks of major forest types in the southeastern Xiaoxinganling Mountains using remote sensing datasets as well as multiyear fieldwork data. We constructed a carbon stock model based on remote sensing data and deep learning, and used the best model to make pixel-by-pixel predictions of carbon stocks to map their spatial distribution time series. This allowed analyzing the spatial and temporal variation characteristics of carbon stocks.

2. Materials and Methods

2.1. Study Areas

The total area of the reserve is 12,133 ha (Figure 1), and it is geographically located at the eastern edge of the Eurasian continent, with temperate continental monsoon climate characteristics [15]. Within the territory of the reserve, the mountains are surrounded by gullies and ravines, forming a more complex mountainous terrain. The average altitude of the whole area is approximately 400 m, and the relative height is between 80 and 300 m. To facilitate forest management and forestry surveys, the forest management department divided the reserve into 444 forest subcompartments according to stand composition and topography (steep slopes, flat land and pit land) of the study area, which made the stand composition and topographic distribution in each forest subcompartments basically the same, and we randomly selected one sample point in each forest subcompartment as the sample point for subsequent research. Because the zoning principle of subcompartments in this study area is basically the same small class with homogeneous forest types, the carbon stock information of each sample point represents the average carbon stock level of the subcompartment.

2.2. Data Acquisition and Treatment

2.2.1. Ground Survey Data

We collected and compiled data from the forest resource planning and design surveys conducted in 1989, 1999, 2009, and 2019, which consisted of comprehensive census data of protected areas within forest subcompartments. The data included tree height (m), diameter at breast height (cm), number of forest trees (n ha⁻¹) and storage volume (m³ ha⁻¹) from 444 forest subcompartments. To obtain the forest biomass of local native tree species more precisely, the biomass per hectare (Mg ha⁻¹) of different forest types was calculated according to the anisotropic growth equations [30] proposed by Dong (Table 1). These growth equations have been extensively experimented and validated, and have high accuracy in the calculation of forest biomass. Based on Yu’s study on the carbon content of forests in Xiaoxinganling, a carbon conversion factor of 0.4412 was applied [31]. Forest biomass and the carbon conversion factor were used to derive forest carbon stocks (Mg C ha⁻¹) for 444 forest subcompartments for 4 years. Due to the excessive number of pixels in the study area, we randomly selected one sample point on each forest subcompartment, which means there are 444 × 4 = 1776 sample points as the target variable. The formula for calculating the forest standing wood carbon stock (CS) is:

$$CS=\left(Y_S+Y_B+Y_F+Y_R\right)\times \rho$$

(1)

Y_S, Y_B, Y_F and Y_R are the biomasses of the stems, branches, foliage and roots of trees, respectively, and $\rho$ (0.4412) is the carbon stock conversion factor.

Table 1. Equations used for estimating the biomass of the different forest stands from the Xiaoxinganling Mountains following Dong 2015 [30]. In these equations, H is the average tree height (m), D is the average diameter at breast height (cm), N is the stand density (n/ha) and V is the storage volume (m³/ha).

Forest Types	Biomass Function
Birch forest	YS = (H/(−0.7161 + 1.7316H))V
	YB = (H/(71.1504 + 2.3594H))V
	YF = (D/(52.1765 + 31.5260D))V
	YR = (D/(−7.7814 + 5.2684D))V
Korean pine forest	YS = (N/(369.6842 + 1.5593N))V
	YB = 0.1807D0.1196H−0.1086V
	YF = 0.2835D−0.2246H−0.2737V
	YR = H/(−7.4947 + 5.4H))V)
Larch forest	YS = 0.2638D0.6084H−0.3052V
	YB = 0.0893D0.5372H−0.5798V
	YF = (D/(−0.320.565 + 69.0763D))V
	YR = (D/(44.2954 + 1.182D))V
Populus forest	YS = 0.2472D0.134H−0.0461V
	YB = (D/(70.9902 + 8.7291D))V
	YF = (D/(−71.0595 + 22.4834D))V
	YR = (D/(−7.72 + 0.1178D))V
Pinus sylvestris forest	YS = 0.3063D0.2638H0.0044V
	YB = 0.0809D0.3271H−0.282V
	YF = (D/(−71.0595 + 22.4834D))V
	YR = (D/(−25.67 + 7.3893D))V
Mixed broad-leaved forest	YS = 0.5193D0.1529H−0.0981V
	YB = (D/(79.2051 + 1.2886D))V
	YF = (D/(−84.3635 + 42.2363D))V
	YR = (H/(−6.079 + 5.4717H))V
Mixed coniferous broad-leaved forest	YS = (D/(−1.3193 + 2.0249D))V
	YB = (D/(58.0685 + 5.6768D))V
	YF = (D/(−507.269 + 78.1762D))V
	YR = (H/(−2.1908 + 5.6433H))V
Mixed coniferous forest	YS = (D/(0.6756 + 2.1D))V
	YB = 0.0631D0.3538H−0.2781V
	YF = (D/(−156.921 + 43.4676D))V
	YR = (D/(14.289 + 4.883D))V

Table 1. Equations used for estimating the biomass of the different forest stands from the Xiaoxinganling Mountains following Dong 2015 [30]. In these equations, H is the average tree height (m), D is the average diameter at breast height (cm), N is the stand density (n/ha) and V is the storage volume (m³/ha).

Forest Types	Biomass Function
Birch forest	YS = (H/(−0.7161 + 1.7316H))V
	YB = (H/(71.1504 + 2.3594H))V
	YF = (D/(52.1765 + 31.5260D))V
	YR = (D/(−7.7814 + 5.2684D))V
Korean pine forest	YS = (N/(369.6842 + 1.5593N))V
	YB = 0.1807D0.1196H−0.1086V
	YF = 0.2835D−0.2246H−0.2737V
	YR = H/(−7.4947 + 5.4H))V)
Larch forest	YS = 0.2638D0.6084H−0.3052V
	YB = 0.0893D0.5372H−0.5798V
	YF = (D/(−0.320.565 + 69.0763D))V
	YR = (D/(44.2954 + 1.182D))V
Populus forest	YS = 0.2472D0.134H−0.0461V
	YB = (D/(70.9902 + 8.7291D))V
	YF = (D/(−71.0595 + 22.4834D))V
	YR = (D/(−7.72 + 0.1178D))V
Pinus sylvestris forest	YS = 0.3063D0.2638H0.0044V
	YB = 0.0809D0.3271H−0.282V
	YF = (D/(−71.0595 + 22.4834D))V
	YR = (D/(−25.67 + 7.3893D))V
Mixed broad-leaved forest	YS = 0.5193D0.1529H−0.0981V
	YB = (D/(79.2051 + 1.2886D))V
	YF = (D/(−84.3635 + 42.2363D))V
	YR = (H/(−6.079 + 5.4717H))V
Mixed coniferous broad-leaved forest	YS = (D/(−1.3193 + 2.0249D))V
	YB = (D/(58.0685 + 5.6768D))V
	YF = (D/(−507.269 + 78.1762D))V
	YR = (H/(−2.1908 + 5.6433H))V
Mixed coniferous forest	YS = (D/(0.6756 + 2.1D))V
	YB = 0.0631D0.3538H−0.2781V
	YF = (D/(−156.921 + 43.4676D))V
	YR = (D/(14.289 + 4.883D))V

2.2.2. Remote Sensing Data

Image data were acquired and processed in the GEE platform [32], and the goal was to collect continuous image data from 1989 to 2021; therefore, we choose data from the Landsat 5 TM (1989–2011), Landsat 7 ETM (2012) and Landsat 8 OLI (2013–2021) surface reflection sensor series images (LANDSAT/LT05/C01/T1_SR, LANDSAT/LE07/C01/T1_SR and LANDSAT/LC08/C01/T1_SR), which have a resolution of 30 m and adequate continuity [33]. The Landsat SR dataset has undergone radiation correction and atmospheric correction. We preprocessed the dataset with cloud removal and fusion processing [12]. The image data time range of 1 June to 31 October each year was chosen to match the ground survey data time and to eliminate the effect of seasonal differences on the data processing results. Finally, time series data with clear pixels were obtained by taking the median value.

The obtained time series images were used to calculate the features needed for the experiment (

Table 2. Original bands, vegetation indices and texture features included in the analysis performed through the GEE platform. The corresponding calculation formula is provided for remote sensing factors when applicable.

Remote Sensing Data	Bands, Indices or Parameters	Definition or Calculation Formula
Original band	BLUE(B1)	Water penetration, distinguishing soil vegetation
	GREEN(B2)	Distinguishing different vegetation
	RED(B3)	Observation of roads, bare soil, and vegetation types
	NIR(B4)	Estimating biomass and distinguishing wet soils
	SWIR1(B5)	Distinguishing the road, bare leaking soil
	SWIR2(B6)	Heat distribution mapping, rock identification
Vegetation Index	Normalized Difference Vegetation Index (NDVI) [35]	(NIR − RED)/(NIR + RED)
	Difference Vegetation Index (DVI) [36]	NIR − RED
	Ratio Vegetation Index (RVI) [37]	NIR/RED
	Atmospheric Ratio Vegetation Index (ARVI) [38]	(NIR − (2 × RED − BLUE))/(NIR + (2 × RED − BLUE))
	Soil-Adjusted Vegetation Index (SAVI) [39]	((NIR − RED)/(NIR + RED + 0.5))(1 + 0.5)
	Weighted Difference Vegetation Index (WDVI) [40]	NIR − 0.5 × RED
	Modified Soil-Adjusted Vegetation Index (MSAVI) [41]	(2NIR + 1 − ((2 × NIR + 1) − 8 × (NIR − RED))0.5)/2
	Modified Normalized Difference Water Index (MNDWI) [42]	(NIR − GREEN)/(NIR + GREEN)
	Land Surface Temperature (LST) [43]	T/(1 + (λT/ρ)lnε); λ is the central wavelength of thermal infrared band, ρ = 1.438 × 10−2 m·k; ε is the surface specific emissivity
	Normalized Difference Building and Soil Index (NDBSI) [44]	(SI + IBI)/2; SI and IBI represent the soil index and building index, respectively
	Wetness (WET) [45]	β1BLUE + β2GREEN + β3RED + β4NIR + β5SWIR1 + β6SWIR2; β1–6 are the different coefficients corresponding to different sensor types
Texture features	Angular Second Moment (ASMB1–6)	The texture metrics are calculated from the grayscale co-occurrence matrix around each pixel in each band. The grayscale co-occurrence matrix is a list of the frequency of occurrence of different combinations of pixel luminance values (grayscale) in an image. It calculates the number of times a pixel with value X is adjacent to a pixel with value Y in a specific direction and distance, and then derives statistics from this table [46].
	Contrast (CONTB1–6)
	Correlation (CORRB1–6)
	Variance (VARB1–6)
	Inverse Difference Moment (IDMB1–6)
	Sum Average (SAVGB1–6)
	Sum Variance (SVARB1–6)
	Sum Entropy (SENTB1–6)
	Entropy (ENTB1–6)
	Difference variance (DVARB1–6)
	Difference entropy (DENTB1–6)
	Information Measure of Corr. 1 (IMCORR1B1–6)
	Information Measure of Corr. 2 (IMCORR2B1–6)
	Max Corr. Coefficient (MAXCORRB1–6)
	Dissimilarity (DISSB1–6)
	Inertia (INERB1–6)
	Cluster Shade (SHADEB1–6)
	Cluster prominence (PROMB1–6)

), including 6 raw bands (blue, green, red, NIR, SWIR1, SWIR2), 11 vegetation indices, and 18 texture features corresponding to the 6 raw bands, for a total of 125 features. Since different features have different magnitudes and units of magnitude, and as different satellite data may produce certain differences that can affect the results of data analysis, we normalized all features to eliminate the effect of magnitude among indicators [34]. Then, we used four years of ground survey data (1989, 1999, 2009 and 2019) as the target variables and chose the stepwise regression method to reduce the dimensionality of the feature variables. Next, we introduced each of the 125 feature variables into the regression equation, conducted significance tests and filtered 25 significant features (GREEN, NIR, SWIR1, ARVI, LST, MNDWI, MSAVI, NDBSI, RVI, SAVI, WDVI, DISS_B1, PROM_B1, VAR_B1, CONT_B2, VAR_B2, DISS_B3, INER_B3, PROM_B3, SVAR_B3, VAR_B3, DVAR_B4, SAVG_B4, DISS_B5, DISS_B6).

3. Methodology

To develop a highly accurate and highly mobile regional scale model for predicting standing forest carbon stocks, we extracted spatial, temporal and spectral information from multitemporal multispectral data [32]. The independent variables were obtained by calculating vegetation indices and texture features, and the significant variables were screened by stepwise regression. In the prediction models, we first selected a multiple linear regression model (MLR) [47], a deep neural network model (DNN) [48] and a convolutional neural network model (CNN) [49] to determine the best model for predicting the standing forest carbon stock. The best predictive model was used to forecast the carbon stocks from 1989 to 2021 pixel by pixel, and specific distribution maps were drawn. A trend analysis was performed on the prediction results to determine the future direction of the standing wood carbon stock from forests. The specific process is shown in Figure 2.

3.1. Models

3.1.1. Multiple Linear Regression (MLR)

First, we selected a multiple linear regression model with parameter visualization for the study, with 25 characteristic variables identified as independent variables and forest carbon stock as the dependent variable. The number of sample points for 4 years was 1776, and 70% of the sample points (n = 1243) were randomly selected for the modeling of carbon stock estimation, while the remaining 30% of the sample data (n = 533) were used for accuracy verification. The mathematical model was:

$${Y=\beta }_0{+\beta }_1{X}_1{+\beta }_2{X}_2{+\beta }_3{X}_3{+\beta }_4{X}_4+\cdots {+\beta }_{\rho }{X}_{\rho }+\epsilon$$

(2)

where Y is the forest carbon stock. X_1,2,3…ρ are the independent variables, where ρ indicates the identity (number) of the variable. β_ρ is the corresponding regression coefficient of each variable, and ε is the model residual.

3.1.2. Deep Neural Networks (DNN)

A deep neural network was used to simulate the distribution of forest carbon stocks in the Liangshui Nature Reserve for 4 years. DNN is a powerful feed-forward machine learning mechanism that provides information about linear and nonlinear associations between a set of independent and dependent variables [50]. DNN is an artificial neural network with multiple layer structures, including input layers, multiple hidden layers and output layers. The hidden layers and built-in hyperparameters of DNN make them superior to other statistical models and classical machine learning models [51]. There is a connective structure between nodes in adjacent layers of the model, but there are no connections between nodes in the same layer or across layers. Different nodes are connected to each other in the form of fully connected layers with a combination of activation functions to build a multilayer complex neural network model. The model structure is shown in Figure 3.

DNN have been widely used in remote sensing in recent years [52]. DNN are trained using forward and backward propagation (BP) algorithms, which continuously update each weight (w) and bias (b) in the network. The backward propagation algorithm, which is the core of the DNN, essentially back-propagates the error between the predicted output value and the true value and continuously modifies the magnitude of w and b in each connection layer of the network model so that the error value is continuously reduced until it meets the accuracy requirement [53]. After the value of the loss function declines to the required threshold through the training of a large amount of data, the value of the parameter at this point is taken as the model parameter, and the regression prediction is performed accordingly. The weight and bias of the BP algorithm are calculated as follows:

$$w*=w-\alpha \frac{\partial J\left(w,b\right)}{\partial w}$$

(3)

$$b*=b-\alpha \frac{\partial J\left(w,b\right)}{\partial b}$$

(4)

where w* and w are the weights after and before the update, respectively; b* and b are the biases after and before the update, respectively; α is the artificially set learning rate; and J(w, b) is the loss function of the model. The optimal parameters w* and b* are obtained through continuous gradient descent iterations.

3.1.3. Convolutional Neural Network (CNN)

A convolutional neural network is a multilayer supervised, learning, fully connected or non-fully connected neural network structure, and the model structure is shown in Figure 4. The structural model of the convolutional neural network includes an input layer, a convolutional layer, a fully connected layer and an output layer [54]. The convolutional layer is a unique structure of the convolutional neural network that contains a set of kernels, or parameters that are learned throughout the training of a pooling layer that contains the data. The input signal and filter are convolved in the convolutional layer, and the implicit features of the input signal are iteratively extracted as the input of the next layer until the mapping relationship in the convolutional neural network converges to the predefined input–output mapping relationship or reaches the termination condition [55]. The number of convolutional layers in the network model is set to 8, and the size of the convolutional kernel is set to 1 × 3. Meanwhile, the ReLU [56] function is chosen as the activation function to avoid gradient disappearance due to the increase in the number of layers and the decline in the learning ability [57]. A max pooling layer is added after every two convolutional layers for downsampling, and only the information with larger weights in the feature mapping is retained to reduce the computational effort [58]. Finally, the information is aggregated to the fully connected layer for information regularization. A total of 25 remote sensing features were input as independent variables, 70% of the sample points were extracted as the training set (n = 1243) and the remaining 30% of the sample data were used as the test set (n = 533) for accuracy verification.

3.2. Evaluation Indicators

3.2.1. Model Evaluation Indicators

In this study, the fitting effectiveness of the MLR, DNN and CNN models was evaluated using three metrics. These metrics include mean absolute error (MAE), root mean square error (RMSE) and coefficient of determination (R²), where smaller MAE and RMSE values and larger R² values indicate a more accurate fit of the models. The formulas of these evaluation indicators are shown below.

$$MAE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left|\left(y_i-\widehat{y_i}\right)\right|$$

(5)

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n(y_i-\widehat{y_i})}{n}}$$

(6)

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

(7)

where n is the number of samples, $y_i$ is the actual observed carbon stock content of sample i, $\widehat{y_i}$ is the predicted carbon stock content of sample i and $\overline{y_i}$ is the actual observed mean value of carbon stock.

3.2.2. Spatial Autocorrelation Analysis

Exploratory spatial data analysis can be used to determine whether there is a correlation between the spatial distribution of standing forest carbon stocks and their neighboring regions and to describe the magnitude of the correlation. The main indicator for calculating spatial autocorrelation is Global Moran’s I [59]. Global Moran’s I is the most commonly used statistic in global correlation analysis and is mainly used to describe the average degree of correlation between a spatial unit and its surrounding area. The formula is [60]:

$$\mathrm{Globalmoran}’\mathrm{s} \mathrm{I}=\frac{N{{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{w}_{ij}\left(x_i-\mu \right)\left(x_j-\mu \right)}{({{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{w}_{ij}){{\displaystyle \sum }}_i{\left(x_i-\mu \right)}^2}$$

(8)

where x_i and x_j are the carbon stock values at i and j, $\mu$ is the mean carbon stock value, w_ij is the spatial weight value between i and j and N is the total number of pixels. The value of Global Moran’s I varies between -1 and 1. When the value of Global Moran’s I is close to −1, it indicates that the negative spatial carbon stock autocorrelation is more significant, while when the value of Global Moran’s I is close to 1, it indicates that the carbon stock clustering effect is more significant. A Global Moran’s I value of 0 indicates no spatial autocorrelation.

3.2.3. Sen’s Slope Estimator

A simple nonparametric procedure developed by Sen (Sen, 1968; Hirsch et al., 1982) was used to estimate possible linear trends in the time series. The formula for calculating Sen’s slope is (10):

$$\beta =\mathrm{median}\left(\frac{X_j-X_i}{j-i}\right) \forall j>i$$

(9)

where β is the slope and the other term of the function represents the median value, where X_j and X_i are the data values at time j and i, respectively.

3.2.4. Mann–Kendall Statistical Test

The Mann–Kendall statistical test (MK) is a nonparametric statistical test that can be applied to the trend analysis of variables that are not normally distributed and are not subject to outlier images. The MK time series test formula is:

$$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}f\left(X_j-X_i\right) \forall j>i$$

(10)

$$f\left(X_j-X_i\right)=\left\{\begin{array}{c}+1 if\left(X_j-X_i\right)>0\\ 0 if\left(X_j-X_i\right)=0\\ -1 if\left(X_j-X_i\right)<0\end{array}\right.$$

(11)

$$Var\left(S\right)=\frac{1}{18}\left[n\left(n-1\right)\left(2n+5\right)-{\displaystyle \sum }_{\mathrm{p}=1}^q{t}_p(t_p-1)\left(2t_p+5\right)\right]$$

(12)

where X_j and X_i are the data values at time j and i (j > i), respectively, n is the length of the time series, p is the number of repetitions, q is the number of unique numbers and t_p is the number of repeats for each repetition. The formula for calculating the standard normal variable Z_MK is:

$$Z_{MK}=\left\{\begin{array}{c}\frac{S-1}{\sqrt{Var\left(S\right)}} if S>0\\ 0 if S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}} if S<0\end{array}\right.$$

(13)

when |Z| ≤ Z_1−α/2, the original hypothesis is accepted, and the trend is not significant; if |Z| > Z_1−α/2, the original hypothesis is rejected, i.e., the trend is considered significant. If the significance level α = 0.05, then the critical value Z_1−α/2 = ±1.96. When the absolute value of Z is greater than 1.96 and 2.58, the trend passes the significance test with 95% and 99% confidence levels, respectively [25]

3.2.5. Hurst exponent

The Hurst index is a key method for describing time series information [28] and it is widely used in meteorology and hydrology. The R/S method is based on the following principle: for a time series X_t (n) (n = 1, 2...), the mean series and the cumulative deviation are calculated as:

$$y\left(\tau \right)=\frac{1}{\tau }{\displaystyle \sum }_{n-1}^{\tau }{X}_t\left(n\right),\tau =1,2, \cdots$$

(14)

$$f\left(n,\tau \right)={\displaystyle \sum }_{n=1}^{\tau }\left(X_t\left(n\right)-y\left(\tau \right)\right),1\le n\le \tau$$

(15)

From this, the polar deviation and standard deviation are calculated:

$$R(\tau )=\underset{1\le t\le \tau }{\max }f(n,\tau )-\underset{1\le t\le \tau }{\min }f\left(n,\tau \right),\tau =1,2, \dots$$

(16)

$$S(\tau )={\left[\frac{1}{\tau }{\displaystyle \sum }_{n-1}^{\tau }{\left(X_t\left(n\right)-y\left(\tau \right)\right)}^2\right]}^{0.5},\tau =1,2,\dots$$

(17)

$$R(\tau )/S(\tau )={\left(c\tau \right)}^H$$

(18)

where $\tau$ is the number of segments that split the time series into segments, n is the length of the time series, X_t(n) is the time series of length n, $y\left(\tau \right)$ is the mean of the time series, $f\left(n,\tau \right)$ is the cumulative deviation of the time series, c is a constant and H is the Hurst exponent, 0 < H < 1. Taking the logarithm of the left and right sides of the above equation yields

$$log(\frac{R(\tau )}{S(\tau )})=1-H log \tau +H log c$$

(19)

The Hurst exponent is denoted by H, and its value ranges from 0 to 1. When 0 < H < 0.5, it indicates that the future carbon stock trend is antiphased to the past trend, referred to as anti-sustainability change, and the closer the H value is to 0, the stronger the anti-sustainability. When 0.5 < H < 1, it indicates that the future carbon stock trend is the same as the past trend, which indicates continuous change in the same direction, and the closer the H value is to 1, the stronger the continuous change. When the value of H is 0.5, it indicates that the carbon stock time series is in an independent state and it constitutes a random series; thus, the future trend and the past trend are not correlated.

4. Results

4.1. Descriptive Statistics

The mean carbon stock value for four years (1989, 1999, 2009, 2019) of field survey data (total 1776 sample points) was 74.44 Mg C ha⁻¹ with a standard deviation of 32.37 Mg C ha⁻¹. The specific annual carbon stock statistics are shown in

Table 3. Ground survey data statistics.

Year	CS-Mean (Mg C ha−1)	CS-Standard Deviation (Mg C ha−1)	CS-Median (Mg C ha−1)
1989	74.42	38.35	71.65
1999	73.62	40.41	67.29
2009	77.25	24.41	78.55
2019	72.46	21.87	70.57

.

4.2. Evaluation of Models

In order to compare the fitting effect of different models, the ratio of training and validation data was chosen as 7:3 (1243:553), and ten combinations were randomly selected for control experiments to finally take the average value to evaluate the models. The evaluation results of the validation data of each model (MLR, DNN, CNN) are shown in

Table 4. The MLR, DNN and CNN model result summary for validation data.

Combination	MLR		DNN		CNN
	MAE	R2	MAE	R2	MAE	R2
1	17.58	0.34	16.87	0.51	16.32	0.54
2	18.16	0.29	17.09	0.48	16.87	0.52
3	18.01	0.35	17.32	0.42	16.34	0.54
4	17.96	0.31	17.29	0.45	17.11	0.49
5	17.60	0.34	16.84	0.51	16.38	0.53
6	17.87	0.32	16.98	0.49	16.90	0.51
7	18.21	0.28	17.13	0.47	16.23	0.57
8	17.54	0.35	17.28	0.46	16.27	0.56
9	17.71	0.32	17.32	0.41	16.97	0.50
10	17.97	0.31	17.19	0.43	16.84	0.52
Mean	17.86	0.32	17.13	0.46	16.62	0.53

. The traditional MLR model fit worse than the other models, while the DNN model fit better. The CNN model has the best fit robustness with carbon storage; the mean MAE is 16.62 Mg C ha⁻¹ and the mean R² is 0.53. Figure 5 shows the correlation between the observed and predicted carbon stocks of the different models, and the MAE decline curve, as well as the loss curve during the training of the CNN model, indicating that the model is not over-fitted. It can be seen that the CNN has a greater advantage in fitting carbon stocks.

4.3. Optimal Model for Prediction

To obtain a more intuitive picture of the spatial and temporal variations in forest carbon stocks from 1989 to 2021, we used the CNN model, which had the highest robustness, to make predictions. All pixel values of the 25 significant variables for each year were substituted into the model to estimate the carbon stock distribution for the whole region in each year, and the results of the specific distribution are shown in Figure 6.

Meanwhile, we smoothed the mean carbon stock fold line (CS-CNN) of the time series using Savitzky–Golay filters [61] to better study the trend of carbon stock changes in the long time series (Figure 7). From the smoothed curve (CS-SG), we observed that the carbon stock tended to be stable from 1989 to 1995 (60–70 Mg C ha⁻¹). However, it rose slightly to the peak from 1995 to 2001 (75–80 Mg C ha⁻¹) and then began to decline. From 2001 onward, the carbon stock started to rise significantly with time and reached a peak in 2007 (90–95 Mg C ha⁻¹). From 2007 to 2012, there was a significant downward trend in carbon stocks. From 2012 to 2021, there was a weak upward trend in carbon stocks (60–70 Mg C ha⁻¹). We applied a trend fit (CS-Trend) to reflect the overall change in carbon stock in the time series, which showed that the overall change in carbon stock in the time series had increased.

4.4. Spatial Autocorrelation Analysis of Carbon Stocks

To investigate whether there are other conditions for the spatial clustering phenomenon in the distribution of standing forest carbon stock in the study area, Moran’s I was introduced for analysis. From the predicted CS-CNN results, 2000 pixels (30 m × 30 m) were randomly extracted for spatial autocorrelation analysis. Because the change in forest carbon stock is a slow process, we selected CNN prediction results every 5 years (1990, 1995, 2000, 2005, 2010, 2015 and 2020) to study the autocorrelation analysis of forest carbon stock, which can avoid the chance error caused by the analysis of carbon stock in a single year. The Global Moran’s I of CS-CNN in 1990, 1995, 2000, 2005, 2010, 2015 and 2020 were 0.12, 0.44, 0.20, 0.14, 0.25, 0.29 and 0.18, respectively (Figure 8), all of which passed the test at the 1% significance level, and the results indicated that the spatial autocorrelation of carbon stock was significant at the 99% confidence level. The Moran scatterplot is divided into four quadrants, which indicate four local spatial correlations. Points in quadrant 1 are high-value regions and are surrounded by high-value regions; points in quadrant 2 are low-value regions and are surrounded by high-value regions; the point in quadrant 3 is a low-value region and is surrounded by low-value regions; the points in quadrant 4 are high-value regions and are surrounded by low-value regions. The Moran’s I of carbon stock in the study area is generally greater than 0, which shows a positive spatial correlation and indicates that carbon stock has different degrees of clustering effects.

4.5. Trends in the Spatial Transformation of Carbon Stocks

Sen’s slope and the MK test were used to detect the spatial variation trends of the carbon stock, and the results were classified into eight categories: extremely significant rise, significant rise, slightly significant rise, no significant rise, extremely significant decline, significant decline, slightly significant decline and no significant decline. The specific distribution and percentages of the change trends are shown in Figure 9 and summarized as follows: 40.9% of the carbon stock showed no significant rise, mainly concentrated in the central part of the study area; 31.7% of the carbon stock showed no significant decreasing trends, mainly concentrated in the south and north of the study area; 18.6% showed extremely significant increasing, significant increasing and slightly significant increasing trends, concentrated in the north, west and east of the study area; 8.8% showed extremely significant decreasing, significant decreasing and slightly significant decreasing trends, mainly observed in the south and north-central regions of the study area. The total increasing trend accounted for 59.5% of all carbon stock, and the total decreasing trend accounted for 40.5%, indicating that the carbon stock in the study area has shown an overall increasing trend over the past 33 years, but some areas have exhibited a decreasing trend.

4.6. Future Trends in Forest Carbon Stocks

To better explore the future trends of carbon stocks in the study area, the Hurst exponent was calculated for the CNN prediction results from 1989 to 2021. This was then combined with the Sen and MK tests to determine the relationship between future evolution and past trends in the changes to vegetation cover in Heilongjiang Province (Figure 10). The Hurst exponent ranged from 0 to 0.80, with a mean value of 0.41, indicating that in the carbon stock study area, the overall change was characterized by reverse persistence since the Hurst exponent range between 0 and 0.5 is set as reverse persistence, and the exponent range between 0.5 and 1 is set as same direction persistence. However, Hurst exponents closer to 0.5 indicate weaker correlations between future changes and past trends for either reverse persistence or continuous change. The spatial carbon stock trends (significant increase, increase, significant decrease, decrease) were superimposed over the sustainability exponent H of these two types of change, and eight development trends were derived, including four types of benign development directions and four types of malignant development directions. Thus, this work demonstrated the change characteristics and future development trends of standing forest carbon stock in the study area.

The different future trends of forest carbon stocks as a percentage are shown in

Table 5. The percentage of different future change trends in CS.

Change Directions	Future Trends	Percentage
Continuous decline	Positive significant decline	0.01
	Positive decline	0.06
Rise in the past but declining trend in the future	Reverse significant decline	0.04
	Reverse decline	0.29
Decline in the past but rising trend in the future	Reverse significant rise	0.13
	Reverse rise	0.41
Continuous rise	Positive significant rise	0.01
	Positive rise	0.05

. The future development trend of the study area as a whole is likely to develop in a benign direction (60% of the future upward trend), and 6% of the regions are continuously rising (1% and 5% of the region show positive significant upward and positive upward trends, respectively). These few regions are distributed in the northwestern and southern parts of the study area. In the study area, 54% of the regions (the largest share of the trend) declined in the past but are on an upward trend now and are likely to remain so in the future. Of these regions, 13% show an anti-significant upward trend and 41% show anti-upward trends. These locations are mainly distributed in the central and western parts of the study region. The regions with a continuous decline accounted for 7% of the study area (positive significant decline and positive declining trend accounted for 1% and 6%, respectively), and were mainly distributed in the central and northwestern parts of the study region. The regions that showed a past rising trend but future declining trend accounted for 33% of the area (reverse significant decline and reverse declining trend accounted for 4% and 29%, respectively), and were concentrated in the north and south parts of the study region.

5. Discussion

Forest carbon stock is a component of the Earth’s carbon pool and an important indicator of ecosystem stability [62]. As a typical representative of the natural boreal secondary forest, the changes observed in the standing forest wood carbon stock in the study area reveal, to a certain extent, the variability of standing forest wood carbon stock in the whole Xiaoxinganling region. Therefore, understanding and predicting the spatial and temporal distribution of standing forest wood carbon stocks is essential for optimizing forest resource management and improving the carbon sequestration potential in this region, as well as mitigating adverse effects of climate change and related impacts.

In this study, a robust CNN model was trained using publicly available and easily accessible Landsat series satellite imagery and existing ground data. Then, the trained model was used to predict the carbon stock distribution from 1989 to 2021. Initially, we tried the MLR model with visible model parameters [63], but the carbon stock of forest stands and the related remote sensing factors did not have a completely linear relationship, and the performance of all indicators was average. Then, we tried deep learning algorithms that could better fit the nonlinear relationship and found that the CNN was the most suitable algorithm for carbon stock prediction, with high portability and better performance for high-dimensional input variables [64]. Because the carbon stock at a location is influenced by environmental conditions and surrounding locations, steep slopes can also affect the distribution of carbon stock [65]. CNN introduces spectral information and texture features, which can reflect the feature of the ground, into the model as independent variables, and also solves the problem of spatial heterogeneity to some extent. However, the CNN also has some limitations, such as underestimating the maximum or minimum carbon stock content due to the limited training set, which is also a drawback of global prediction. In subsequent research, local regression combined with deep learning can be used to make more accurate predictions of carbon stock content.

To observe the dynamic changes in continuous carbon stock time series, we choose the image datasets from the Landsat series of satellites [66], which have the longest operation and are continuously maintained. However, there are still some gaps between the remotely sensed data and the observations made directly in the field due to atmospheric attenuation, the low signal-to-noise ratio of the sensor itself, annual revisits and low resolution [67]. In particular, as the scale of the study area is extended from local to regional, the degradation of image quality will lead to greater errors in the analysis. In addition, defects in the remote sensing data itself can interfere with the experiment. In this study, remote sensing images from some years (2001, 2012) showed slight streaking due to the failure of the on-board scan line corrector of the Landsat satellite. The appearance of streaks in the original bands led to the same situation in the remote sensing factors extracted from them, which led to similar streaks in the CS-CNN predicted forest carbon stock mapping (Figure 6). The stripes makes the mapping look less than perfect, but the mean values of forest carbon stock projections are in the normal range and have little effect on the trend analysis (Figure 7). The other years of images do not have this problem, and the satellite data from 1989, 1999, 2009 and 2019 used to train the model are also normal, which means that the robustness of the model itself is not affected in any way. Subsequent observations and studies can include higher-resolution satellite hyperspectral images to predict or verify the changes in carbon stock patterns and multisource remote sensing data can be combined to achieve better results [68].

The results of this study show that the overall trend of forest carbon stocks in the study area is increasing. Especially between 1989 and 2007, the increasing trend was clear. Under the reasonable management of the protected area, the growth of forest stands reached the mature phase, and the stumpage carbon stock was close to the peak. Afterwards, most of the stands entered the overmature forest stage, the resistance of trees weakened and the carbon stock began to decline due to natural disasters before finally stabilizing. As the resistance of the overmature forest weakens, rapid climate change and the influence of human factors led to a decreasing trend in forest carbon stocks. Ultimately, forest carbon stock stabilized due to the self-regulation ability of the forest ecosystem.

The study found that forest carbon stock does not continually increase or decrease by following a single trend but increases or decreases periodically, which is consistent with natural relationships and reveals the feasibility and practicality of the study to a certain extent. However, we need to conduct more in-depth research to quantitatively explore the drivers of forest carbon stock changes, which includes anthropogenic factors and climate change. These considerations will help to more accurately predict forest carbon stock trends, which will in turn facilitate more precise forest carbon stock distribution mapping.

6. Conclusions

Based on the GEE platform and Landsat TM/ETM/OLI remote sensing data, we studied the long-term spatial and temporal variations in the carbon stocks of forest stands using remote sensing and deep learning methods and explored the development trend of carbon stocks in the study area in depth. The MLR model (RMSE = 24.00, R² = 0.34) had a general week fitting effect. Therefore, we chose a deep learning algorithm to predict the carbon stocks in forest stands, and we found that this model (RMSE = 17.45, R² = 0.64) had the best fit. The CNN with the best robustness was used to make pixel-by-pixel predictions of the carbon stock. It was found that the carbon stock was mostly stable from 1989 to 1995; after 1995, the carbon stock started to rise significantly and then began to decline and stabilize after reaching its peak in 2007. There was a weak upward carbon stock trend between 2012 and 2021. Then, we analyzed the trend of the carbon stock of forest stumpage in recent years by combining Sen and MK statistical tests. The total upward trend was 59.5%, and the total downward trend was 40.5%. Generally, the forest stumpage carbon stock in the study area has developed positively in recent years. Combined with the Hurst index, we analyzed the future development trend of the standing wood carbon stocks of the forests in the study area and found that the overall development trend for 60% of the study area is likely to be benign in the future. Subsequently, we will use more source data and more comprehensive methods to further study the characteristics of the changes in the entire northern ecosystem.