Temporal and spatial dynamics of Net Primary Productivity and prediction of wetland carbon sequestration potential on the Tibetan Plateau

BP neural network prediction model

In 1986, a research team led by Rumelhart and McClelland proposed the BP neural network (Shavlik & Towell, 1989). This structure is a multilayer feedforward network trained using an error backpropagation algorithm, representing one of the most widely used neural network models currently (Chai & Kun, 2021). BP neural networks possess the advantages of approximating arbitrary nonlinear mapping relations, enhanced generalization ability, improved fault tolerance, and simple implementation (Wang, 2015).

Core formulas

(1) Hidden layer output: (2.1)${\displaystyle h_j=f\left(\sum _{i=1}^n{w}_{ij}{x}_i+b_j\right).}$ (2) Output layer output: (2.2)${\displaystyle y_k=g\left(\sum _{j=1}^m{w}_{jk}{h}_j+b_k\right).}$ Where:

x_i: Input variables; h_j: Output of the j-th hidden layer neuron; y_k: Model predicted value.

w_ij: Weight between the i-th input and j-th hidden neuron; w_jk: Weight between the j-th hidden neuron and k-th output neuron; b_j: Bias of the j-th hidden neuron; b_k: Bias of the k-th output neuron.

f(⋅): Hidden layer activation function (commonly Sigmoid function: $f\left(x\right)=\frac{1}{1+e^{-x}}$; g(⋅): Output layer activation function (linear function for regression tasks).

Model accuracy verification

(1) Correlation coefficient (R)

For two variables, x andy, the sample correlation coefficient R is defined as follows: (2.3)${\displaystyle R=\frac{\sum _{i=1}^n\left(x_i-\stackrel{-}{x}\right)\left(y_i-\stackrel{-}{y}\right)}{\sqrt{\sum _{i=1}^n{\left(x_i-\stackrel{-}{x}\right)}^2}\sqrt{\sum _{i=1}^n{\left(y_i-\stackrel{-}{y}\right)}^2}}}$ where n is the sample size; x_i and y_i are the ith observation of variables x and y, respectively; and $\stackrel{-}{x}$ and $\stackrel{-}{y}$ are the sample means of x and y, respectively.

(2) Model robustness validation (Monte Carlo simulation)

A Monte Carlo perturbation simulation was conducted to assess the robustness of the neural network model in the presence of input uncertainties. The procedure is summarized as follows:

(a) Input perturbation

The original input matrix was represented as follows: (2.4)${\displaystyle p=\left\{p_{ij}\right\},i=1,\dots ,m,j=1,\dots ,n}$ where m is the number of features, and n is the number of samples.

At the k- th simulation, each input element can be perturbed as follows: (2.5)${\displaystyle p_{ij}^{\left(k\right)}=p_{ij}\cdot \left(1+\rho \cdot \left(2r_{ij}-1\right)\right)}$ where ρ = 0.02 is the perturbation ratio (±2%), and r_ij ∼ U(0,1) is a random number uniformly distributed in [0,1].

(b) Model prediction

The perturbed input $p^{\left(k\right)}$ was normalized and then fed into the trained neural network model f(): (2.6)${\displaystyle y^{\left(k\right)}=f\left(p^{\left(k\right)};\theta \right),k=1,\dots ,N}$ where θ is the fixed network parameter. The number of simulations was set to N = 300.

(c) Statistical analysis

The outputs $\left\{y^{\left(1\right)},\dots ,y^{\left(N\right)}\right\}$ were aggregated for statistical evaluation:

Mean prediction: (2.7)${\displaystyle \stackrel{-}{y}=\frac{1}{N}\sum _{k=1}^N{y}^{\left(k\right)}.}$ Standard deviation (uncertainty): (2.8)${\displaystyle {\sigma }_y=\sqrt{\frac{\sum _{k=1}^N{\left(y^{\left(k\right)}-\stackrel{-}{y}\right)}^2}{N-1}}.}$ Mean squared error (MSE) distribution: (2.9)${\displaystyle MS{E}^{\left(k\right)}={\frac{\sum _{j=1}^n\left(y_j^{\left(k\right)}-{\hat{y}}_j\right)}{n}}^2}$ By calculating the mean and standard deviation of $MS{E}^{\left(k\right)}$, the robustness of the model under input perturbations can be quantitatively assessed. Moreover, plotting the prediction mean alongside the confidence band (mean ±1 standard deviation) provides a visual representation of the variability in model outputs.

(3) Normalized root mean square error (NRMSE) (2.10)${\displaystyle NRMSE=\frac{RMSE}{y_{\max }-y_{\min }}\times 100\text{\%}}$ (2.11)${\displaystyle RMSE=\sqrt{\frac{\sum _{i=1}^n{\left(\widehat{y_i}-y_i\right)}^2}{n}}}$ where n is the number of samples; y_i is the ith observed (true) value; ${\hat{y}}_i$ is the ith predicted value; y_max is the maximum of the observed values; and y_min is the minimum of the observed values.

Kriging interpolation

In this study, the ordinary Kriging interpolation method was employed, while the spherical function was selected as the variational function for the interpolation calculation. In contrast to inverse distance weighting and spline interpolation, Kriging interpolation is based on the variational function model. This approach effectively utilizes the spatial autocorrelation information of the data and enables the optimal unbiased estimation of unknown points (Tang & Yang, 2021).

Assuming that the target interpolation point is p₀, number of sample points is known to be p₁, p₂, ⋅⋅⋅, p_m, and observation corresponding to the sample point is f(p₁), f(p₂), ⋅⋅⋅, f(p_m), the estimate $\tilde{f}\left(p_0\right)$ at p₀ obtained through Kriging interpolation is calculated as follows (Tang & Yang, 2021): (2.12)${\displaystyle \tilde{f}\left(p_0\right)=\sum _{k=1}^m{w}_{k}f\left(p_k\right)}$ where w_k is the weight coefficient and satisfies ${\sum }_{k=1}^m{w}_k=1$.

CASA modeling

The CASA model was developed by a team of researchers from the Carnegie Institution, Stanford University, and NASA’s Ames Research Center (Potter et al., 1993). Precipitation, temperature, radiation data, and Normalized Vegetation Index (NDVI) data for the same period were obtained for the span of 2025–2030 for this scientific work. The NDVI data were predicted using a BP neural network. Based on the CASA model, the net primary productivity (NPP) of the target study area was simulated so as to analyze the dynamics of the ecosystem productivity in the area and reveal its evolution in a long time series.

In the CASA model, the calculation of net primary productivity (NPP) of vegetation follows this approach: the value of net primary productivity of vegetation for a given image element x at moment t is equal to the product of the photosynthetically active radiation absorbed by that image element x at moment t and the actual light energy utilization of that image element x at moment t. It is expressed in a mathematical formula as: (2.13)${\displaystyle NPP\left(x,t\right)=APAR\left(x,t\right)\times \varepsilon \left(x,t\right)}$ (2.14)${\displaystyle APAR\left(x,t\right)=SOL\left(x,t\right)\times FPAR\left(x,t\right)\times 0.5}$ (2.15)${\displaystyle FPAR\left(x,t\right)=0.5\times \left(FPA{R}_{NDVI\left(x,t\right)}-FPA{R}_{SRi,min}\right)}$ (2.16)${\displaystyle FPA{R}_{NDVI\left(x,t\right)}=\left(NDVI\left(x,t\right)-NDV{I}_{i,min}\right)\times \frac{FPA{R}_{\max }-FPA{R}_{\min }}{NDV{I}_{i,max}-NDV{I}_{i,min}}+FPA{R}_{\min }}$

(2.17)${\displaystyle FPA{R}_{SR\left(x,t\right)}=\left(SR\left(x,t\right)-S{R}_{i,min}\right)\times \frac{FPA{R}_{\max }-FPA{R}_{\min }}{S{R}_{i,max}-S{R}_{i,min}}+FPA{R}_{\min }}$ (2.18)${\displaystyle SR\left(x,t\right)=\frac{1+NDVI\left(x,t\right)}{1-NDVI\left(x,t\right)}.}$ In this equation, NPP(x, t) represents the net primary productivity of the vegetation of pixel x at the moment of t, which has the unit of gC⋅m⁻²⋅a⁻¹; APAR(x, t) is the photosynthetically active radiation absorbed by the pixel at the moment, which has the unit of MJ⋅m⁻²⋅a⁻¹; FPAR_NDVI(x,t) is the proportion of incident photosynthetically active radiation absorbed by the vegetation layer calculated from the NDVI; SOL(x, t) represents the total solar radiation at spatial location x and time t; FPAR(x, t) generally refers to the proportion of photosynthetically active radiation of vegetation; SR(x, t) is the ratio vegetation index of image x in month t; FPAR_SR(x,t) is the proportion of incident photosynthetically active radiation absorbed by the vegetation layer calculated from the ratio vegetation index; NDVI(x, t) is the NDVI of image x in month t; NDVI_i,min and NDVI_i,max are the minimum and maximum values of NDVI for the ith planting type, respectively; SR_i,min and SR_i,max are the minimum and maximum values of SR for the ith planting type, which are calculated from NDVI_i,min and NDVI_i,max for the corresponding ith planting type, respectively; FPAR_max and FPAR_min are the percentages of incident photosynthetically active radiation absorbed by the vegetation layer, respectively. FPAR_max and FPAR_min are the maximum and minimum proportions of incident photosynthetically active radiation absorbed by the vegetation layer, which are 0.95 and 0.001, respectively; and ɛ(x, t) while represents the actual light energy utilization of the pixel x at the moment of t, which has the unit of gC⋅Mj⁻¹. (2.19)${\displaystyle C\left(x,t\right)=NPP\times B.}$ C represents the carbon sequestration potential of the wetland, NPP is the net primary productivity of the wetland vegetation, and B is the carbon conversion factor of the wetland, which takes the value of 0.44 g (Shen et al., 2022).

Transfer matrix

The land use transfer matrix is typically employed to characterize the conversion relationships among various land use types between the initial and final time points within a specific time interval (Zhang et al., 2017).

Core formulas (2.20)${\displaystyle T=\left[\begin{array}{cccc}{T}_{11}\hfill & T_{12}\hfill & \cdots \hfill & T_{1n}\hfill \\ T_{21}\hfill & T_{22}\hfill & \cdots \hfill & T_{2n}\hfill \\ ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ T_{n1}\hfill & T_{n2}\hfill & \cdots \hfill & T_{nn}\hfill \end{array}\right]}$ Core element:T_ij represents the area/quantity of the ith land use type at t₁ transferred to the jth land use type at t₂;

Diagonal elementT_ii: The area/quantity of land use type that remains unchanged from t₁ to t₂;

Off-diagonal elementT_ij:The area/quantity of land use type that undergoes conversion from t₁ to t₂.

Evaluation of model accuracy

In this study, monthly regression analysis was conducted on the NDVI of the Tibetan Plateau from 2001 to 2010. A total of 367,511 monthly data points were utilized over the 10-year period. The analysis yielded an average R-value of 0.7 across the 12 months, an NRMSE of 10.26%, and an MSE of 0.04. The range of the R value was [−1, 1]. An R-value closer to 1 indicates a strong positive correlation, while one nearing −1 signifies a strong negative correlation. The MSE and NRMSE range was [0, +∞); values approaching 0 indicate greater consistency between the model’s predicted and true values. Therefore, the model exhibited good accuracy. Concurrently, Monte Carlo simulations were performed, yielding the following values: Original MSE: 0.0089; Mean MSE after perturbation: 0.0128; Std of MSE: 0.0000. In the context of Monte Carlo simulation perturbation, model error increased; however, the fluctuation remained minimal, thereby indicating good stability (Fig. 4).

Analysis of NDVI spatial and temporal dynamics

The spatial distribution of monthly NDVI on the Tibetan Plateau of China from 2025 to 2030 was characterized by high values in the southeast and low values in the northwest. Notable differences in NDVI between these regions were particularly significant in January, May, and September (Figs. 5–8).

The overall spatial distribution of NDVI on the Tibetan Plateau over the next six years showed relatively high NDVI values in June–August and relatively low values in February, October, and November.

The folded plot of monthly mean NDVI changes on the Tibetan Plateau from 2025 to 2030 indicated that the highest monthly mean NDVI values occurred in August for 2025, 2026, and 2030 and June for 2027–2029. The monthly average NDVI minimums for 2025–2028 occurred in December, while those for 2029 and 2030 were recorded in January (Fig. 9). The highest monthly average NDVI value over the six years was 0.331 in June 2028, while the lowest was 0.064 in December 2027.