A Factor Analysis Backpropagation Neural Network Model for Vegetation Net Primary Productivity Time Series Estimation in Western Sichuan

Abstract

Vegetation net primary productivity (VNPP) is the main factor in ecosystem carbon sink function and regulation of environmental processes. However, NPP data products have data missing in some areas, which affects the availability and overall accuracy level of data. Therefore, we adopted the Factor Analysis Backpropagation neural network model (FA-BP model) to acquire a high-accuracy and high-reliability NPP result without missing or empty areas by using a series of easily accessible datasets, such as meteorological data and remote sensing data. We selected the western Sichuan region as the study area and carried out a VNPP time series estimation from 2000 to 2016. Comparative simulations also verify the accuracy of the time series estimation results: The Pearson correlation r of VNPP prediction results ranged from 0.807 to 0.917, the mean absolute error ranged from 29.1 to 38.9, the root mean square error was between 37.3 and 51.8, and the mean relative error varies from 0.10 to 0.14. Further analysis shows that the spatial pattern of estimated VNPP during the past 17 years in western Sichuan shows a decreasing trend from southeast to northwest. Besides, the VNPP time series is generally on an upward trend in this period. The increasing and decreasing areas of VNPP values in the study area accounted for 81.42% and 18.58%, respectively. Moreover, we find that temperature dominates the change of VNPP in the whole western Sichuan region. The data processing method and experimental results presented in this paper can provide a reference for accurately acquiring VNPP and related studies on natural resources and climate change.

1. Introduction

VNPP is an essential aspect of ecosystem services, reflecting the productive capacity of vegetation communities in terrestrial ecosystems, characterizing ecosystem quality, determining ecosystem carbon source/sink functions, and evaluating ecological regulation processes [1]. In many ways, surveying measurements is the most reliable way to estimate VNPP with high accuracy. However, the application is limited, the portability is poor, and long-term time series data are difficult to obtain [2]. Most commonly used NPP estimation methods are vegetation estimation models, including statistical models, light energy utilization models, and process models [3]. These models can estimate NPP from their respective advantages, but they all have shortcomings. The statistical model is simple in principle, but the modeling theory is lacking, and the portability is poor. The light energy utilization model requires fewer parameters, is simple and easy to operate, and can be estimated on a large scale. But the shortcoming is the lack of research on the internal mechanism of vegetation [4]. The process model has high simulation accuracy, but the disadvantage is that the necessary parameters are complex and challenging [5]. Therefore, NPP estimation has certain uncertainties no matter which model is invoked, so introducing new estimation methods is unavoidable.

In the seasonal estimation of long time series, understanding the composition of climate change factors, making full use of temporal details, and acquiring high-quality satellite data products play an important role [6]. With the advantages of a large detection range, fast data acquisition, low ground constraints, and a large amount of information, remote sensing images have been one of the most effective tools for monitoring the Earth’s surface. They are widely used in land cover classification, urbanization, and crop monitoring [7,8,9]. The spatial and temporal resolution of remote sensing data are balanced, but high spatial resolution products often lack concentrated and reliable time series. In addition, cloud cover, atmospheric conditions, and sensor errors can also lead to impaired temporal resolution [10]. Therefore, there is inevitably missing data in the estimation process, which leads to impaired estimation accuracy. It greatly limits the application of remote sensing images and increases the difficulty of temporal analysis. To reconstruct the missing data of optical images and to improve the data availability and multi-temporal analysis, we have to smooth the time series data to fill in the missing values. The traditional decomposition of time series data is to build a combination of algorithms by grouping components together in the addition or multiplication phase [11]. Then, scholars established a new method named X11 Decomposition to slow down the seasonal component changes. It suggests a trend cycle for all-time series data but cannot handle the outliers effectively [12]. Yet now, wavelet variation is widely used in the trend significant analysis of time series by combining the decomposition of time series with the description of time series variation. Ghaderpour performed a consistency analysis of temperature time series from six geodesics very long baseline interferometry by the spectral and wavelet analyses techniques, showing how temperature affects the measurement results within the annual scale [13]; Verbesselt uses the breaks for the additive season and trend (BFAST) to decompose the normalized difference vegetation index (NDVI) time series into three components, including trend, seasonality, and residual components. He takes full account of confounding factors and identifies the physical changes independent of the physical indicators [6]; Andronis used the Kalman smoothing filter to process NDVI time series data and compared with Cook’s distance and Hat matrix regression analysis to remove the influence of missing values, and to obtain reliable results of forest weather changes [14]; Zhou proposed a deep-learning-based method to reconstruct the Landsat images of time series, and the error was effectively reduced, which proved the effectiveness of his method [15]; In the estimation of VNPP, it is common to use MOD17A3 NPP data products to estimate VNPP and analyze the spatial and temporal patterns of climate change [16]. It is the annual-scale NPP data calculated by the BIOME-BGG model which considers a variety of parameters. However, NPP data products have missing data in some respects, which affects the data’s availability and overall accuracy level. Therefore, removing the effects of missing values and obtaining complete and reliable VNPP estimates is important, and can provide reliable help for subsequent analysis.

Machine learning algorithms have recently been a hot topic in studying geoscientific problems in recent years. Feng Xia used an artificial neural network model to predict the primary productivity of Daya Bay [17]; Haijuan used a neural network to estimate the NPP of the Guanzhong-Tianshui Economic Zone [18]; Chuanhua estimated the NPP of grassland in the permafrost zone of the Qinghai-Tibet Plateau based on machine learning [19]; Jing used the FA-BP neural network to effectively predict the environmental carrying capacity of coastal areas [20]; and Wenying used the FA-BP neural network to establish an inverse model of wheat leaf powdery mildew, which improved the estimation accuracy and had good compatibility [21]. The FA-BP neural network is widely used with simple principles and can be used to deal with common problems that require a combination of various factors and are difficult to specify. It has good generalization ability, fast convergence, and high accuracy approximation ability and is widely used in various regression prediction and classification problems [22].

This paper firstly combined NDVI, digital elevation model(DEM), and meteorological data with the FA-BP neural network to construct an NPP estimation method to make up for the lack of high-precision data in some areas and to estimate the VNPP accurately and quickly, followed by validation of method validity and accuracy through comparison and simulation of data products. Lastly, we quantitatively analyzed VNPP and spatio-temporal patterns in western Sichuan from 2000 to 2016 and explored the reaction of climate factors with them to provide a theoretical basis for evaluating and constructing ecological and environmental protection in the region.

2. Materials and Methods

2.1. Materials

2.1.1. Study Area

The study area (Figure 1) is located in western Sichuan Province, China, and includes a total of 26 counties and cities in Ganzi Tibetan Autonomous Prefecture and Aba Tibetan and Qiang Autonomous Prefecture, with a total area of 200,928.687 km² (28°03′N~34°18′N, 97°23′E~104°10′E). The region is part of the Qinghai-Tibet Plateau and the Hengduan Mountains in China and is a sensitive and initiating zone for climate change. In addition, it has various climate types, significant geographical differences in vegetation types, and complex topography decreasing from west to east, where there is a high mountain plateau alpine triple-high climate zone, with elevations mostly above 3000 m. Additionally, its noticeable climate change, rich and diverse ecosystems, and vital ecological status make it suitable for VNPP research.

2.1.2. Data SOURCE and Preprocessing

NDVI is a radiometric value reflecting the relative abundance and activity of living green vegetation. It is often used to characterize the physiological condition of vegetation, green biomass, and vegetation productivity in the study area [23] through image stitching, reprojection, and format conversion of MODIS NPP and NDVI data using MODIS-specific processing tools(modis reprojection tool (MRT)). The outliers were removed, and a reasonable range of data was selected and multiplied by a scale factor to convert the grayscale values to actual values to obtain usable MODIS NPP and NDVI data. To eliminate extraneous environmental influences, the NDVI data were processed by the maximum value composite (MVC) [24] to obtain the NDVI time series raster datasets from 2000 to 2016. The meteorological data was a dataset of monthly values of terrestrial climate information in China, which underwent strict quality control and met the accuracy and completeness requirements of the study. The monthly values were assembled into annual scale data, and the meteorological data were interpolated using Anusplin professional interpolation software and statistical product service solutions software (SPSS). DEM data were introduced as covariates to eliminate the influence of geographic location and features on spatial interpolation and improve interpolation accuracy (

Table 1. Data List.

Data	Name	Time Resolution	Spatial Resolution	Source
NPP	MOD17A3	1a	500 m	NASA
NDVI	MOD13A1	16d	500 m	NASA
Meteorological	Meteorological	\	\	CMSDC
DEM	DEM	\	90	RAEDCP

). The log file and standard deviation image analysis errors showed that the interpolation results of the meteorological station points meet the experimental requirements.

The spatial resolution of data was sampled uniformly to 1km. To ensure that the projected area remains equal to the field, it was projected as Alberts projection by setting Central Meridian, Standard_Parallel_1, Standard_Parallel_2 to 105, 25, 47, respectively, and the geographic coordinate system to WGS 1984 [25]. Then, it was extracted by mask to obtain the NPP, NDVI, DEM, and meteorological raster dataset in western Sichuan.

2.2. Methods

2.2.1. FA-BP Neural Network

The FA-BP neural network is constructed by performing factor analysis on the characteristic parameter data and then training the prediction through the BP neural network (Figure 2). Factor analysis is a statistical technique that studies how to extract common factors from a population of variables. It can distill variables with specific correlations into a smaller number of factors and use them to represent the original variables [26,27]. It can realize the mining of influence factors behind multivariate and dimensionality reduction before data modeling and has better interpretation than principal component analysis. This method can simplify the network structure, accelerate the network’s convergence speed, and help improve the network’s generalization ability with the following equation.

$$X_i=a_{i1}{F}_1+a_{i2}{F}_2+\cdots +a_{im}{F}_m+{\epsilon }_i,\left(i=1,2,\cdots ,p\right)$$

(1)

where $F_1,F_2,\cdots ,F_m$ are called common factors and ${\epsilon }_i$ is a specific factor of $X_i$. The matrix can represent the model as:

$$X=AF+\epsilon$$

(2)

$$X=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_P\end{array}\right],A=\left[\begin{array}{ccc}{a}_{11}& \cdots & a_{1m}\\ ⋮& \ddots & ⋮\\ a_{p1}& \cdots & a_{pm}\end{array}\right],F=\left[\begin{array}{c}{F}_1\\ ⋮\\ F_m\end{array}\right],\epsilon =\left[\begin{array}{c}{\epsilon }_1\\ ⋮\\ {\epsilon }_p\end{array}\right]$$

(3)

and meet$m\le p$; $Cov\left(F,\epsilon \right)$ = 0; $D_F=D\left(F\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \ddots & 0\\ 0& 0& 1\end{array}\right]$; and $D_{\epsilon }=D\left(\epsilon \right)=$ $\left[\begin{array}{ccc}{\sigma }_1^2& 0& 0\\ 0& \ddots & 0\\ 0& 0& {\sigma }_p^2\end{array}\right]$.

2.2.2. NPP Estimation and Analysis Algorithm

The principle of the number of factors extraction is when the cumulative contribution of the characteristic roots reaches 85% or more; or, if the last of the eigenvalues corresponding to the Mth factor are greater than 1, then the first M factors can be used instead of the overall variable.

BP neural network consists of two parts, “BP” means backpropagation, and the neural network represents a complex computational network. As one of the most widely utilized network models, BP comprises the input layer, the hidden layer, the output layer, and the weights between each layer. The input layer is responsible for receiving external signals, and the hidden layer performs internal information processing, mainly for information transformation. The output layer receives the information transformed by the network function and further processes it to get the output result. The weights are responsible for connecting the layers of information.

The algorithm mainly allows the information to enter the network from the input layer up to the output layer and then operates continuously back and forth through error backpropagation. In forward propagation, the input signal passes through the implicit layer after nonlinear transformation to obtain the output signal:

$$S_j={{\displaystyle \sum }}_{i=0}^{m-1}{w}_{ij}{x}_i+b_j$$

(4)

$$x_j=f\left(S_j\right)$$

(5)

where$w_{ij}$ is the weight between nodes and nodes, $b_j$ is the threshold value of the node j, $x_i$ is the output value,$f$ is the activation function, generally selected as an S-shaped function or linear function.

If there is a discrepancy with the actual value, then the error backpropagation process is transferred to adjust the weights and thresholds from the output to the input direction:

$$E\left(w,b\right)=\frac{1}{2}{{\displaystyle \sum }}_{j=0}^{n-1}{\left(d_j-y_j\right)}^2$$

(6)

where $d_j$ are all the results of the output layer.

Repeated learning training until the desired result is approximated. According to the gradient descent method, the weights and thresholds between the implicit and output layers are adjusted.

$$w_{ij}=w_{ij}-{\eta }_1·\frac{\partial E\left(w,b\right)}{\partial w_{ij}}$$

(7)

$$b_j=b_j-{\eta }_2·\frac{\partial E\left(w,b\right)}{\partial b_j}$$

(8)

The adjustment of the weights and thresholds between the input and implicit layers are as follows.

$$w_{ki}=w_{kj}-{\eta }_1·\frac{\partial E\left(w,b\right)}{\partial w_{ki}}$$

(9)

$$b_i=b_i-{\eta }_2·\frac{\partial E\left(w,b\right)}{\partial b_i}$$

(10)

where $w_{ki}$ is the weight between nodes and nodes,$b_i$ is the threshold value of the node j, and${\eta }_1$ and${\eta }_2$ represent study efficiency.

The trained neural network can process similar sample information on its own to output information with minimal error after nonlinear transformation [28,29,30].

To build an inversion model for VNPP time series estimation, control points were evenly arranged in the high-quality inversion area (NPP quality control(NPP_QC) ≤ 32) within the study area. Among them, 70% of the data were randomly selected as training samples and the remaining 30% as test samples. As is shown in Figure 3, from the DEM, meteorology, and NDVI data, a total of 10 characteristic parameters were selected to drive the model in this paper, including 10 characteristic parameters such as elevation, slope, slope direction, and so on. This paper performed factor analysis on the sample data retained the principal components for dimensionality reduction, and reduced the 10 feature parameters to 5 to form an input data set jointly with NPP data products. Then, a neural network model was constructed with five input layers, four hidden layers, and one output layer. The best network model was obtained and used for VNPP inversion by adjusting the number of training times and learning rate, setting the transfer function, constantly changing the weights and thresholds to reduce errors, and repeatedly testing the reliability of the output results. The inversion results were then compared with NPP data MOD17A3 and the results of other similar studies for validation.

The spatio-temporal pattern of VNPP is obtained by trend and Rescaled Range analysis. Then, the response of climate factors to VNPP is received by related analysis.

(1)

Trend Analysis Method

To quantitatively describe the trend and intensity of VNPP, the trend analysis method was selected to conduct a linear regression of the state-by-image calculation of the trend propensity rate Slope for VNPP from 2000–2016 [31]. Slope > 0 means that the NPP has increased over the 17 years, whereas the opposite is a decrease. In addition, the larger the absolute value, the more significant the change. A statistical test of correlation coefficient was used for the significance trend test. The changing trend is considered at a significant or highly significant level if the confidence level of the correlation coefficient is 95% or 99% or more. Accordingly, the trends of change can be classified into five categories: highly significant increase (0 < p < 0.01, slope > 0), significant increase (0.01 < p < 0.05, slope > 0), insignificant change (p > 0.05), significant decrease (0.01 < p < 0.05, slope < 0), and highly significant decrease (0 < p < 0.01, slope < 0).

(2)

Rescaled Range Analysis Method

The fractal characteristics and long-term memory process of the VNPP time series are quantitatively described, and the Hurst index calculated on this basis can better measure the statistical correlation of the time series [32]. The correspondence between Hurst index and time trend is as follows: when 0 < H ≤ 0.5, the time series is anti-persistent or inverse state persistent, that is, the future direction is opposite to the past, and the closer to 0, the stronger the anti-persistence; when H = 0.5, the time series is a standard random wandering state, and it can be considered that the present trend will not have an impact on the future, and there are no long-term correlation characteristics; when 0.5 ≤ H < 1, it means that the time series is an ongoing series and there is state persistence, i.e., the future trend is the same as the past and the closer to 1, the stronger the persistence.

(3)

Related Analysis

The influence of single climatic factors (temperature, precipitation, and solar radiation) on VNPP was investigated. The correlation coefficients between each climatic factor and VNPP were calculated using image element-based correlation spatial analysis [33]. A correlation coefficient > 0 indicates a positive correlation and <0 shows a negative correlation, and the larger the absolute value, the stronger the correlation. In addition, the t-test was used to test the significance of the correlation coefficients.

3. Results and Analysis

3.1. Validation of NPP Estimation Results

The KMO and Bartlett’s sphericity tests were conducted for the fitness of the factor analysis results. The test results (

Table 2. Results of the factor analysis adaptation test.

Time	KMO	AC
2000	0.665	16,248.5
2001	0.767	14,025.1
2002	0.718	22,273.4
2003	0.750	15,745.3
2004	0.651	21,999.7
2005	0.687	13,404.2
2006	0.738	17,596.8
2007	0.769	18,793.5
2008	0.743	18,161.1
2009	0.685	10,111.9
2010	0.686	15,886.2
2011	0.700	13,158.3
2012	0.699	10,501.4
2013	0.722	17,819.9
2014	0.684	14,333.8
2015	0.681	18,221.0
2016	0.799	13,791.2

) showed that the KMO is 0.651~0.799 are more significant than 0.6, the approximate chi-square is 13,404.2~21,999.7, the degrees of freedom are all 45, and the significance levels are all 0.000. Therefore, the original hypothesis is invalid, and factor analysis can be used with better results.

The correlation coefficients between the NPP and BP neural network for vegetation in western Sichuan from 2000 to 2016 were calculated (

Table 3. Correlation coefficients of the training set, test set, validation set, and BP neural network.

Time	Training	Validation	Test	All
2000	0.86512	0.82594	0.89347	0.86244
2001	0.87986	0.89885	0.91339	0.88679
2002	0.90095	0.88478	0.90217	0.89872
2003	0.88292	0.84220	0.88838	0.87851
2004	0.88652	0.88888	0.86180	0.88086
2005	0.88062	0.82577	0.85725	0.87168
2006	0.87618	0.85087	0.84339	0.86681
2007	0.89335	0.86147	0.90641	0.89092
2008	0.91130	0.85679	0.90599	0.90252
2009	0.83066	0.81766	0.82209	0.82703
2010	0.87184	0.81120	0.85792	0.86246
2011	0.86306	0.88205	0.86997	0.86674
2012	0.85335	0.86516	0.80543	0.84798
2013	0.91872	0.85889	0.91189	0.90904
2014	0.88753	0.83188	0.84608	0.87030
2015	0.88271	0.87930	0.87950	0.88197
2016	0.89166	0.86096	0.86535	0.88385

). Among them, the training set was used in the modeling process; the validation set was used to adjust the model to ensure the “training” effect has a good generalization ability to ensure the success of the “test.” The test set was used for the evaluation of the final model to assess the model. The test set was used to evaluate the final model and assess the model’s generalization ability and only works once. In the training sample, the percentage of the training set, test set, and validation set = 70:15:15. The correlation coefficients of the BP neural network with each data set of training samples and MODIS NPP simulated products are 0.85335~0.91872, 0.81120~0.89885, 0.84339~0.91339, and 0.82703~0.90904 and are more significant than 0.8, indicating that the BP neural network fit is better.

The model was tested for simulation accuracy using the test sample, the independent variables of the test sample were input, and the inverse simulation of VNPP was performed using the BP neural network model established above. The results were analyzed for errors (

Table 4. The error of simulated values with MOD17A3 NPP product.

Time	Fitting Equation	r	p-Value	MAE	RMSE	MRE
2000	y = 0.9734x − 0.5042	0.882	0.008	34.4	44.6	0.12
2001	y = 1.0365x + 0.6028	0.824	0.009	32.3	41.3	0.11
2002	y = 0.9801x − 2.9982	0.831	0.007	32.5	40.7	0.10
2003	y = 0.9890x + 1.6989	0.887	0.007	33.1	42.1	0.12
2004	y = 0.9784x + 1.0849	0.856	0.005	29.1	37.3	0.10
2005	y = 0.9729x + 3.6431	0.827	0.006	33.5	41.6	0.10
2006	y = 1.0867x − 35.362	0.917	0.008	38.9	51.8	0.14
2007	y = 1.0046x − 3.8315	0.815	0.007	35.3	45.2	0.13
2008	y = 1.0347x − 23.968	0.803	0.006	35.5	44.5	0.12
2009	y = 0.9928x − 0.8811	0.896	0.006	35.8	44.9	0.12
2010	y = 1.0092x − 4.5111	0.842	0.007	33.6	43.8	0.11
2011	y = 0.9869x + 2.0913	0.857	0.005	34.3	44.4	0.12
2012	y = 0.9757x + 4.4076	0.913	0.005	33.4	41.5	0.11
2013	y = 0.8937x + 1.5209	0.824	0.007	32.1	40.3	0.11
2014	y = 0.9935x − 5.2277	0.836	0.009	34.1	43.2	0.12
2015	y = 1.0509x − 4.4948	0.807	0.008	35.2	46.2	0.12
2016	y = 1.0174x − 4.0790	0.876	0.008	38.4	53.8	0.14

). The Pearson correlation coefficient r of the VNPP prediction results from 2000–2016 ranged from 0.807 to 0.917, and the results had a good correlation. All p-values were less than 0.01, meeting the statistical requirements (Figure 4). The mean absolute error ranged from 29.1 to 38.9, the root mean square error ranged from 37.3 to 51.8, and the mean relative error ranged from 0.10 to 0.14, all within a reasonable range. The slope of the fitted equation between MODIS NPP products and predicted values ranged from 0.8937 to 1.0867, close to 1. The prediction results were evenly distributed and did not show any overall bias, indicating that the model has good applicability.

Comparative analysis and reference to the previous research results were made to revalidate the results of this study. The VNPP estimated in this study is close to Xinqian [34] using MOD17A3 NPP data products for VNPP estimation in southwest China. The overall results are slightly higher than those Siyao [35] estimated for the whole province of Sichuan based on the Carnegie-Ames-Sanford approach model (CASA). The Qinghai-Tibet Plateau region within Sichuan province in this study is roughly the same as the study by Chuanhua [36]. The tested estimation results are as follows (Figure 5).

3.2. Spatial Pattern Characteristics of NPP

Influenced by climate change, human activities, etc., the spatial pattern of the annual average VNPP in western Sichuan from 2000 to 2016 varies significantly (Figure 6). The overall spatial pattern decreases from southeast to northwest, ranging from 42.7 to 705.1 g C·m⁻²·a⁻¹, with a mean value of 375.16 g C·m⁻²·a⁻¹. The overall level of VNPP in western Sichuan is moderate, with below-average areas widely distributed. These areas are mainly concentrated in the west part of Ganzi Tibetan Autonomous Prefecture and the northwestern part of Aba Tibetan Autonomous Prefecture, accounting for 56.8% of the region’s total area. The climate mainly belongs to the Qinghai-Tibetan plateau climate. The terrain is stepped down from west to east, the topography is a hilly plateau and high plain, and the vegetation type is mainly an alpine meadow [37]. Shiqu County is located at the intersection of Tibet, Sichuan, and Qinghai provinces, with an altitude above 4000 m. The terrain is a primarily hilly plateau with significant diurnal temperature differences, long sunshine, and low temperatures [38]. Therefore, VNPP is the lowest in the study area, with an average value of only 187 g C·m⁻²·a⁻¹. The Dege-Sedar-Ganzi-Xinlong-Litang-Daocheng area is consistent with the direction of the Yalong River flow. The VNPP is comparable to the mean value, between 218 and 523 g C·m⁻²·a⁻¹, and is extremely rich in hydropower resources, good foliage, and dense forests. Because of altitude, the eastern part of Rangtang and Aba County areas are influenced, with significant differences in temperature and precipitation.

The vegetation is dominated by mixed forests of spruce-fir and subalpine evergreen coniferous forests. The high-value VNPP areas are mainly located in the southwestern part of Aba, the southeastern part of Ganzi, and a small amount in the southwest part of Ganzi Autonomous Prefecture and the western part of Baiyu County, all near the Sichuan Basin or Sichuan Mountain area, with hydropower resources and minerals. The vegetation cover types are alpine meadows and alpine swamp meadows. Among them, the site from Jiuzhaigou to Kangding is in the transition zone from the western Sichuan plateau and mountains to the Sichuan Basin, and the highest VNPP ranges from 350.6 to 705.13 g C·m⁻²·a⁻¹. The average value is up to 512.2 g C·m⁻²·a⁻¹. The climate is influenced by circulation and complex topography with significant vertical climates, forming a monsoon climate. The area is rich in biological resources, and the vegetation is mostly grassland, scrub, and forest. The VNPP in Ruo’ergai-Hongyuan-Maerkang-Jinchuan-Daofu-Yajiang is between 315.4~643.6 g C·m⁻²·a⁻¹ and the temperature and precipitation are significantly affected by altitude. This area is one of the critical areas of biodiversity in China, along with the advance of the western Sichuan plateau area towards the Sichuan mountain region. The climate gradually changes from a highland cold-temperate humid climate to a Tibetan plateau sub-humid climate. In addition, local areas such as Hongyuan, Kangding, and Xiangcheng counties significantly differ in VNPP NPP from their surroundings, indicating significant regional differences in VNPP.

3.3. Characteristics of NPP Changes

The time series trend of annual mean VNPP from 2000–2016 varies significantly (Figure 7), and the fluctuation range of annual mean VNPP in western Sichuan during the study period is 338.0777~409.2726 g C·m⁻²·a⁻¹, with the lowest value being 2000, and the highest value being 2006. The annual average VNPP in the whole region shows a fluctuating upward trend with an increased rate of 1.528, passing the significance test (p < 0.01). The spatial variation characteristics of the VNPP vary significantly (Figure 8), and the overall trend shows spatial heterogeneity. The VNPP of the study area overall increases during this period, and the rooms with decreasing and increasing trends are 18.58% and 81.42%, respectively. Only 8.14% passed the significance test, and 91.86% of the areas showed insignificant changes. Among them, the significant increase (p < 0.05, Slope > 0) accounted for 7.34% of the site, mainly in the Hongyuan and Ruoerge grassland areas, central Ganzi Autonomous Prefecture, and around Yajiang County, which are rich in natural resources and biodiversity. In recent years, human activities have become increasingly influential, resulting in increased vegetation cover, improved ecological environment, and a clear trend of increased NPP [39]. The significant decrease (p < 0.05, Slope < 0) accounts for only 0.8% of the area, mainly in some Jiuzhaigou, Maerkang, and Xiangcheng counties. The vegetation growth in these areas has been degraded, and the VNPP trend has decreased significantly due to geological hazards and climate change.

In addition to the influence of climatic conditions on VNPP, human activities are also an unavoidable factor. In recent years, with the expansion of the scope of human activities, infrastructure construction, farmland reclamation, bridge, and road construction have caused the decline of VNPP. However, at the same time, the structure of ecological environment protection, reforestation, and the prohibition of overgrazing are also essential in the recovery of VNPP [40,41].

The persistence of VNPP trends was further analyzed (Figure 9). The Hurst index was calculated by the R/S analysis method, with a balance of 0.13~0.86 and a mean value of 0.43. The areas below and above the mean value accounted for 48.83% and 51.17%, respectively, indicating that the western Sichuan region has a weak persistence. Nearly half of the trend will not continue its current development trend. Combining Slope and Hurst index analysis, the development favorable area (Slope > 0, hurst > mean or Slope < 0, hurst < mean) should insist on ecological environment construction work to maintain good stability of the local ecosystem. The development unfavorable areas (Slope < 0, hurst > mean or Slope > 0, hurst < mean) should take environmental protection measures in time to prevent the decline of the vegetation growth trend.

3.4. Impact Factor Analysis

Climate factors play a dominant role in influencing VNPP [42]. This paper analyzed the correlation between temperature, precipitation, solar radiation, and VNPP on an image-by-image basis using year as the time scale (Figure 10). The correlation coefficients of the three climate factors and VNPP are −0.62 to 0.86, −0.86 to 0.67, and −0.63 to 0.82, respectively. There is some spatial variation in the response of VNPP to climate factors.

The VNPP is mainly positively correlated with temperature, with 74.61% of the positively correlated areas, of which 18.70% are significant areas, primarily concentrated in the western part of the study area, where high altitude, low temperature, and insufficient heat resulted in limited vegetation growth. Therefore, there is an apparent positive correlation between the two. VNPP and precipitation are mainly negatively correlated, with 81.54% of the negatively correlated areas, including 16.16% of the significant places found in all western Sichuan regions. Western Sichuan has a semi-humid, semi-arid climate. It is at a high altitude, making it difficult for warm and humid air currents to cross into the plateau, resulting in relatively little rainfall but concentrated in the summer and autumn months [43]. Because of the complex topography of western Sichuan, the over-concentrated precipitation forms continuous heavy rain and leads to the frequent occurrence of various geological floods. It threatens the survival environment of vegetation and thus inhibits the growth of VNPP, so the two have a stronger negative correlation. VNPP and solar radiation are mainly positively correlated, with 82.67% of the positively correlated areas, of which 14.63% are significant. The high altitude of the western Sichuan area, sufficient sunshine hours, high intensity of sunshine, and high atmospheric transparency do not have enough weakening effect on radiation. Therefore, the amount of solar radiation reaching the ground is more [44]. Adequate radiation prolongs the growth time of plants and is more conducive to plant photosynthesis and increased output, so the two have a noticeable positive correlation [45].

In terms of significance, the temperature is dominant in influencing the VNPP in western Sichuan.

4. Discussion

This paper takes the lead in introducing the FA-BP neural network model to build a novel VNPP estimation method, which is available to compensate for the loss of accuracy caused by the lack of NPP data. Moreover, we proved that it is sufficient to support VNPP estimation by combining easily accessible data sets (i.e., NDVI, DEM, and meteorological data) and that the acquired results are accurate and reliable.

Most researchers mainly used light energy and process models for VNPP estimation, but both have different shortcomings. Since light energy use models lack mechanistic studies within plants, previous VNPP estimation methods have certain limitations. For example, among the many parameters of the CASA model for estimating VNPP inputs, photosynthetically active radiation (PAR) and light energy utilization are determinants. They have a powerful influence on the estimation accuracy, and the parameters are influenced not only by vegetation type but also by the uniformity of vegetation cover, elevation, slope, etc. [46,47]. Process models are complex in structure, and model parameters are not easily accessible. For example, there are uncertainties in the parameter information of the boreal ecosystem productivity simulator model (BEPS model), including meteorology, soils, and ground cover, and also require flux observations of different ecosystems for validation [48,49]. The principle of the FA-BP neural network [50] is simple, does not need to consider complex ecological processes, is easy to program and implement, and is faster and more practical than the previous two. In addition, in comparison with other similar studies, the FA-BP neural network can effectively reduce the number of output parameters, thus increasing the training rate of the model and weakening the contradiction between training rate and accuracy, which is faster and more practical than other machine learning algorithms.

The novel method proposed in this paper is available to obtain the long-time VNPP series. The single requirement for input parameters is uniform spatio-temporal resolution, regardless of the spatio-temporal influence of specific topography. Therefore, it has good robustness and is capable of different application scenarios. Furthermore, the final output is a high-precision estimation result without missing data area, which supports the complete spatio-temporal pattern and climate factor response research in a large study area.

5. Conclusions

This paper introduced an FA-BP neural network model for VNPP time series estimation for solving the missing data problem in NPP data products. Combining various easily accessible data products, NDVI, DEM, and meteorological data, a VNPP time series estimation and high spatio-temporal resolution results without blank or empty areas were acquired.

To further understand the ecological environment changes in the western Sichuan area over the past decade, this paper used a trained model to acquire the VNPP changes and its spatial-temporal pattern analysis in western Sichuan from 2000 to 2016. The VNPP time series analysis showed that the average VNPP in western Sichuan ranges from 42.7 to 705.1 g C·m⁻²·a⁻¹, with a multi-year average of 375.16 g C·m⁻²·a⁻¹, decreasing from southeast to northwest. Furthermore, the VNPP values in the study area showed an overall fluctuating upward trend, with 81.42% and 18.58% of the study area experiencing increases and decreases, respectively. In addition, the response of climate factors to VNPP was explored. The analysis showed that the role of climatic factors depends on whether the hydrothermal conditions suitable for plant growth are met, and the dominant factor of VNPP changes in western Sichuan is temperature.

The FA-BP neural network model successfully estimated VNPP results with high reliability and accuracy, which helps to solve missing data problems for existing NPP data. In addition, the proposed method has good robustness, which is suitable for different application scenarios. The estimated high-precision NPP data can reasonably evaluate the effectiveness of ecological construction work in the region in recent years. Moreover, it provides an essential theoretical basis for exploiting natural resources and climate change studies and gives a reference for the area to make sustainable development policies. However, there is a contradiction between the learning rate and stability of FA-BP neural networks. And it lacks an effective method to determine the number of neurons. Further research could consider speeding up training while avoiding entering local minima to ensure stability for improving the performance of the method.