Net Primary Productivity Estimation of Terrestrial Ecosystems in China with Regard to Saturation Effects and Its Spatiotemporal Evolutionary Impact Factors

: The net primary productivity (NPP) of vegetation holds a pivotal character for the global carbon balance as a key parameter for characterizing terrestrial ecological processes. The most commonly used indices for estimating vegetation NPP, for instance, the normalized difference vegetation index (NDVI), often suffer from saturation issues that can compromise the accuracy of NPP estimation. This research utilizes a new vegetation index based on the radial basis function (RBF) to estimate vegetation NPP in Chinese terrestrial ecosystems over the past two decades (2001–2020) and investigates the spatiotemporal variation characteristics of NPP and the driving mechanisms. The results indicate that the kernel vegetation index ( k NDVI) can effectively alleviate the saturation problem and signiﬁcantly improve the accuracy of NPP estimation compared to NDVI. Over the past two decades, the NPP of Chinese terrestrial vegetation ranged from 64.13 to 79.72 g C/m 2 , with a mean value of 72.75 g C/m 2 , showing a ﬂuctuating upward trend. Changes in the NPP of terrestrial ecosystems in China are mainly affected by precipitation. The dominant factors inﬂuencing NPP changes varied over time and had different impacts. For instance, in the period of 2001–2005 the climate had a positive effect on NPP changes, with the dominant factors being evaporation and precipitation. However, in the period of 2010–2015 the dominant climate factors shifted to evaporation and temperature, and their effect on NPP changes became negative. The outcomes of this research aim to serve as a foundation for carbon cycle research and ecosystem environment construction in China.


Introduction
Net primary productivity (NPP) is the remaining portion of the total organic matter produced by photosynthesis per unit area of green vegetation per unit time, excluding its autotrophic respiration [1]. It straightforwardly indicates the productivity and quality of an ecosystem in the context of natural ambient conditions, and is one of the significant elements in deciding the carbon source/sink of an ecosystem [2,3]. In the current context of global climate change, the accurate and timely estimation of the NPP in China's land ecosystems is critical to implementing the "Double Carbon Strategy" and the advancement of ecological civilization in the country [4].
The integration of remotely sensed data and ecological models has emerged as the primary method for large-scale vegetation NPP estimation owing to its broad coverage, low leading to lower uncertainty in the propagation of kernelized vegetation indices compared to non-kernelized ones.
Therefore, in this study, a new vegetation index is used to address the saturation phenomenon in estimating the vegetation NPP of Chinese terrestrial ecosystems over the past two decades, then the Theil-Sen median and Mann-Kendall trend are used to explore the trends of vegetation NPP in Chinese terrestrial ecosystems. Second-order partial correlation, multiple correlation analysis, the differential equation method, and the PLS-SEM method are used to analyze the role of each influencing factor on the change in vegetation NPP in Chinese terrestrial ecosystems. This study can provide a valuable scientific reference for ecological environmental conservation, carbon cycle, and resource allocation research in China.

Research Region
China is located in the southeastern part of Eurasia, along the western coast of the Pacific Ocean, and is characterized by a strong thermal contrast between land and sea [26,27] (Figure 1). China possesses the world's most typical monsoon climate zone, with high temperatures and heavy precipitation in summer and cold and poor precipitation in winter, when the high-temperature period coincides with the rainy period [28]. China is home to a rich diversity of vegetation types, ranging from cold boreal coniferous forests to warm subtropical broadleaf evergreen forests and tropical rainforests. The main vegetation types include broadleaf evergreen forests, deciduous broadleaf forests, coniferous forests, grasslands, and deserts [29]. These unique characteristics make China an important area for ecological research and natural resource management.

Material and Processing
(1) Remote-sensing data For the research, the remote sensing data we required included vegetation-type data and NDVI data. Vegetation-type data images with a spatial resolution of 500 m ( Figure 2) were retrieved from h ps://lpdaac.usgs.gov (accessed on 10 August 2022) and were maintained by the NASA EOSDIS Land Processes Distributed Activity Archive Center (LP DAAC) at the USGS Earth Resources Observation and Science (EROS) Center (h ps://doi.org/10.5067/MODIS/MCD12Q1.006, accessed on 10 August 2022). The vegetation-type data were then reclassified into thirteen categories according to the needs of the

Material and Processing
(1) Remote-sensing data For the research, the remote sensing data we required included vegetation-type data and NDVI data. Vegetation-type data images with a spatial resolution of 500 m ( Figure 2) were retrieved from https://lpdaac.usgs.gov (accessed on 10 August 2022) and were maintained by the NASA EOSDIS Land Processes Distributed Activity Archive Center (LP DAAC) at the USGS Earth Resources Observation and Science (EROS) Center (https://doi.org/ 10.5067/MODIS/MCD12Q1.006, accessed on 10 August 2022). The vegetation-type data were then reclassified into thirteen categories according to the needs of the study ( Table 1). The spatial distribution dataset of China's monthly 1 km vegetation index (NDVI) was based on the ten-day 1 km vegetation index data of SPOT/VEGETATION PROBA-V 1 km PRODUCTS (http://www.vito-eodata.be, accessed on 6 June 2022) using the maximum value synthesis method.  (Table 1). The spatial distribution dataset of China's monthly 1 km vegetation index (NDVI) was based on the ten-day 1 km vegetation index data of SPOT/VEGETATION PROBA-V 1 km PRODUCTS (h p://www.vito-eodata.be, accessed on 6 June 2022) using the maximum value synthesis method.  (2) Meteorological data This research required the following meteorological data: monthly average temperature, monthly total precipitation, monthly total solar radiation, and potential evapotranspiration data. The temperature, precipitation, and evaporation data were obtained from National Earth System Science Data Center and National Science and Technology Infrastructure of China (h p://www.geodata.cn, accessed on 2 September 2022), with a spatial resolution of 1 km. The solar radiation data were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) (h ps://doi.org/10.24381/cds.68d2bb30, accessed on 30 August 2022), with a spatial resolution of 1 km.
(3) NPP validation data The NPP validation data required for the study included both station data and NPP product data. The data of eight carbon flux observation stations were downloaded from the China National Ecological Science Data Centre Resource Sharing Service Platform  (2) Meteorological data This research required the following meteorological data: monthly average temperature, monthly total precipitation, monthly total solar radiation, and potential evapotranspiration data. The temperature, precipitation, and evaporation data were obtained from National Earth System Science Data Center and National Science and Technology Infrastructure of China (http://www.geodata.cn, accessed on 2 September 2022), with a spatial resolution of 1 km. The solar radiation data were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) (https://doi.org/10.24381/cds.68d2bb30, accessed on 30 August 2022), with a spatial resolution of 1 km.
(3) NPP validation data The NPP validation data required for the study included both station data and NPP product data. The data of eight carbon flux observation stations were downloaded from the China National Ecological Science Data Centre Resource Sharing Service Platform (http://www.cnern.org.cn, accessed on 30 October 2022); the specifics of the eight sites are listed in Table 2, with the shared flux data comprising the net ecosystem exchange (NEE) and ecosystem respiration (R e ). To validate and evaluate the reliability of the NPP simulation results, this study converted the flux tower observations from January to December 2010 to NPP (GPP = R e − NEE, NPP EC = α × GPP) [30]. The NPP simulation results of the sites were extracted based on the latitude and longitude information, and a linear correlation analysis was used for NPP validation. The NPP product data used for  (4) Additional data Additional data considered for the study comprised administrative boundary data, soil data, socio-economic data, and topographical data. Data on administrative divisions were taken from the National Centre for Basic Geographical Information (http://ngcc.sbsm.gov. cn/, accessed on 15 June 2022). GDP and population data were taken from the Resource and Environmental Science and Data Centre (https://www.resdc.cn/, accessed on 19 September 2022). Soil texture data were acquired from the Harmonized world soil data base v1.2 of the Geographical FAO database of the Chinese Academy of Sciences. The 30 m resolution DEM data were downloaded from the Geospatial Data Cloud (http://www.gscloud.cn/, accessed on 15 June 2022) and slope data were extracted on the basis of these data (Table 3).

Methods
The technical route for this study encompassed the four main aspects illustrated in Figure 3: (1) Data Acquisition, which involves obtaining remote-sensing data, meteorological data, validation data, and other relevant data sources; (2) Construction of the kernel

Methods
The technical route for this study encompassed the four main aspects illustrated in Figure 3: (1) Data Acquisition, which involves obtaining remote-sensing data, meteorological data, validation data, and other relevant data sources; (2) Construction of the kernel function by referencing the RBF algorithm in machine learning to optimize the normalized vegetation index NDVI; (3) Estimation and Accuracy Validation of the net primary productivity (NPP) utilizing the CASA model to estimate vegetation NPP in Chinese terrestrial ecosystems in the past two decades and validating the results using the validation data; and (4) Analysis of the spatiotemporal evolution of NPP, including trend analysis and examination of the driving mechanisms.

Construction of kNDVI
With reference to the RBF machine learning algorithm [25], this study used linear algebra to perform operations and derive a nonlinear algorithm from the linear algorithm to ensure that the normalized vegetation index NDVI could be optimized.
This study refers to [25] and uses the RBF kernel function k(a,b) = exp(−(a − b) 2 /((2σ 2 )) in all cases, where the σ parameter controls the concept of the distance between the IR and NIR bands. This kernel function is calculated as follows: The kNDVI is calculated by taking the length scale parameter σ to the average distance between the NIR and IR bands with σ = 0.5 (n + r). Thus, the equation for kNDVI can be further simplified as provided below. kNDVI = tanh NDVI 2 (2)

NPP Estimation Model
Based on the modified CASA model from [31], the parameters required for the models were processed accordingly to finally estimate the vegetation NPP with a temporal resolution of one month and a spatial resolution of 1 km. The model estimation equation is where NPP(x,t) is the NPP of pixel x in month t (g C/m 2 ·month 1 ), APAR(x,t) is the photosynthetically active radiation absorbed by pixel x in month t (g C/m 2 ·month 1 ), and ε(x,t) is the actual light energy utilization of pixel x in month t (g C/MJ).

NPP Trend Analysis
The combination of the Theil-Sen median trend analysis with the Mann-Kendall trend test approach is a valuable method to identify trends in long-term sequence data [32]. This method has become increasingly popular for analyzing long-term vegetation data. For this research, we utilized the Theil-Sen median to examine the trends in vegetation NPP over the past two decades in Chinese terrestrial ecosystems.
Theil-Sen median trend analysis is a reliable non-parametric statistical approach to trend analysis that can reduce the impact of outliers [33]. Its formula is where median is the median function and NPP j and NPP i are the NPP in year j and year I, respectively. When β > 0, it means that NPP shows an increasing trend, while the shows a decreasing trend.
The Mann-Kendall test is a non-parametric test that neither demands a sample to obey a specific distribution nor is influenced by a few outliers [34]. During the test, the vegetation NPP results for the period from 2001-2020 were constructed as a set of time series with per-pixel metavalues to determine significance differences: where Z is the statistic of the NPP series, S refers to the test statistic, var(S) is the variance of the S statistic, n denotes the time series length, sgn denotes a symbolic function, and NPP i and NPP j are the time series datasets. In the double-lateral trend test, Z > 0 shows an uptrend, Z < 0 shows a downtrend, and Z = 0 shows unchanged. Thus, provided the level of importance α, if |Z| > u 1−α/2 this implies the fact that the NPP time series shows significant changes at the α level. In this paper, we determined the importance of change trends in vegetation NPP time series at the 0.05 confidence level.

Response of NPP to Climatic Factors
Vegetation NPP is usually influenced by climate elements such as temperature, precipitation, and radiation. Therefore, in this study, temperature, precipitation, and radiation were used as climatological influences on NPP, as well as correlation and significance levels among vegetation NPP, with temperature and precipitation investigated on pixel-by-pixel bases using partial and multiple correlation analyses.
The second-order partial correlation coefficient allows the correlation between the remaining two variables to be analyzed by excluding the interference of two of the four variables [35]. It is calculated as shown below: where r ab,cd denotes the partial correlation coefficients for variables a and b with the constant variables c and d, and r ab,c , r ad,c , and r bd,c denote the partial correlation coefficients for variables a and b, a and d, and b and d, respectively. A t-test was used to test the significance of the statistics: where r 12,34, . . . ,m indicates the bias correlation coefficient, n indicates the sample size, and m indicates the quantity of independent variables. In this study, p ≤ 0.05 was accepted as statistically meaningful. Multiple correlation analysis focuses on the influence of multiple factors on one factor [36]. Assuming that y and z are independent variables and x is the dependent variable, the multiple correlation coefficient between them is calculated as follows: where x is the time series NPP, y represents temperature, and z represents precipitation.
An F-test was conducted as well: where r x,yz is the multiple correlation coefficient, n is the sample size, and k is the independent variable number. In order to discriminate the impacts of climatic factors on vegetation NPP in different areas, the relative response of climatic factors to NPP was classified into five categories based on previous studies (Table 4). Table 4. Relative types of impacts of temperature and precipitation on NPP.
Driving Factor F C P T P P

Strong common influence
F C is the result of the F-test of the multiple correlation coefficient between NPP and temperature and precipitation; P T and P P are the p-values of the partial correlation coefficients between NPP and temperature and precipitation, respectively.

The Role of Climate and Human Activities in NPP Changes
Overall, temperature, precipitation, and solar radiation were used as climate-influencing factors for NPP, and differential equations were used to assess the contributions of climate change and human activities to the vegetation NPP of the terrestrial ecosystems of China for 2001-2020. On each stage, the drivers of vegetation NPP change were discussed using a partial least squares-structural equation (PLS-SEM) approach with a five-year period, integrating four types of influencing factors-meteorological, topographic, soil, and socioeconomic-to reveal how vegetation NPP responded to climatic factors across terrestrial ecosystems of the country.
(1) Differential equation method Differential equation methods can quantitate the response for climate change and human activities to changes in vegetation NPP for a function y = f (x 1 , x 2 , . . . ). The variability of the dependent factor y may be represented as [37]: It then follows that: Alternatively, where d NPP /d t is the trend in the long time series NPP and d NPP /d t , d tem /d t , d pre /d t , and d radi /d t represent the slope of the linear regression of NPP, temperature, precipitation, and solar radiation at time t, respectively. Moreover, C cc = C tem + C pre + C radi , C ha = ε, C tem , C pre , and C radi are the response of climate change, human activity, temperature, precipitation, and radiation to the long-term trend in NPP, respectively, and ε is the systematic error.
Although there are many drivers other than climate change that affect NPP changes, human activities constitute the dominant part. Finally, the individual proportional contributions of climate change (Equation (16)) and human activity (Equation (17)) to NPP trends, p(x), can be estimated as follows: (2) PLS-SEM method The advantages of PLS-SEM are that it fully exploits the information in the data, minimizes the error term, and does not require a particular (high) sample size, model identification problems, or distribution state of the data. Therefore, the problem of covariance between variables can be dealt with effectively. The model is composed of two major components; one is a measuring model characterizing the relationship between latent and explicit variables, and the other is a structuring model characterizing the relationship of latent and explicit variables [38]. Its methodology for estimating the parameterization is separated into two procedures: (1) iteration to obtain the estimates of the latent variables, and (2) linear regression using partial least squares to obtain the estimates of the parameters for structural and measurement models. In this study, using a stratified random sampling method, 50,000 sample points were selected to analyze the effect of each influencing factor on the changes in NPP.
Considering k latent variables and each of k sets of dominant variables which each contain m variables, each dominant variable set may be denoted as It is assumed that the latent variables are all linearly combined with the latent variables and the latent variables are linearly combined with the dominant variables. As each dominant variable is correlated as a distinct latent variable, the equation for measuring the model is The equation of the structural model is as follows: where ξ i is the standardized latent variable, λ ij is the factor loading, β ij is the path coefficient, and σ ij and ε i are error correction terms with mean zero and that are uncorrelated with the predictor.

Accuracy Validation of NPP
This research compared simulated NPP values with data from eight flux sites during the same period. In addition, it compared NPP simulations for forests, shrubs, and grasslands in 2005-2010 with time series data from the Dinghu Mountain, Haibei, and Dangxiong sites. To further evaluate the NPP estimates for various vegetation types, this study included three NPP products (MOD17A3HGF, GLO_PEM, and Chinanpp).
The study results show that the simulated NPP values based on the CASA model are strongly correlated with the measured values at the stations (R 2 > 0.6) ( Figure 4). Interestingly, the NPP estimated with kNDVI (R 2 = 0.74) had a better correlation than that estimated with NDVI (R 2 = 0.65). This finding suggests that the NPP values estimated using kNDVI were more accurate than those obtained using NDVI. Compared with NDVI, the NPP values estimated based on kNDVI were closer to the measured NPP values in both densely vegetated and sparsely vegetated sites ( Figure 5), and more closely matched the phenological cycle measured at the flux stations, making this approach more suitable for estimating terrestrial vegetation NPP on a national scale.     [39] demonstrated that the NPP values estimated by the GLO-PEM model were generally higher than the MOD17A3HGF product NPP values and [40] produced the Chinanpp dataset, with the maximum light energy utilization determined as 0.55, while the maximum light energy utilization estimated by the other two products and in this study was defined according to each vegetation type (Table 1); thus, the NPP values of Chinanpp products were low. The NPP mean estimated by NDVI in this study was high compared with the NPP mean of the MOD17A3HGF product, which was overestimated, while the NPP mean estimated by kNDVI was close to the NPP mean of MOD17A3HGF, and both of them were in better concordance. Generally, the NPP values estimated using kNDVI in this study were reliable and could be used in subsequent analysis. On a monthly scale, the NPP trend line for the study area over the studied period from January to December generally showed an inverted "V" shape ( Figure 8), with the NPP peaking in July each year. The monthly variation features of NPP are similar to the climate characteristics, with higher NPP values in summer and lower NPP values in winter. Because the changes in vegetation NPP are influenced by the control of vegetation growing process and light energy utilization, in turn, light energy utilization is affected by different circumstantial factors, which include temperature, precipitation, and solar radiation. The hot and humid environment in summer is suitable for vegetation growth, while the cold and poor rainfall in winter is not conducive to the accumulation of vegetation organic ma er.   On a monthly scale, the NPP trend line for the study area over the studied period from January to December generally showed an inverted "V" shape ( Figure 8), with the NPP peaking in July each year. The monthly variation features of NPP are similar to the climate characteristics, with higher NPP values in summer and lower NPP values in winter. Because the changes in vegetation NPP are influenced by the control of vegetation growing process and light energy utilization, in turn, light energy utilization is affected by different circumstantial factors, which include temperature, precipitation, and solar radiation. The hot and humid environment in summer is suitable for vegetation growth, while the cold and poor rainfall in winter is not conducive to the accumulation of vegetation organic matter.

Spatial Distribution Characteristics
In the study area, there was a clear pa ern of regional variation in the multi-annual NPP averages, with a general spatial trend showing a decrease from the southeast coast to northwest inland (Figures 9 and 10). High NPP regions were focused in the eastern and southern portions of the mainland, mainly in southwestern Hainan Province, Fujian Province, and most of Taiwan Province. These regions have low latitudes, sufficient heat, and abundant precipitation, which are conducive to the accumulation of organic compounds in vegetation. Low NPP regions were clustered in the western and northern parts of the

Spatial Distribution Characteristics
In the study area, there was a clear pattern of regional variation in the multi-annual NPP averages, with a general spatial trend showing a decrease from the southeast coast to northwest inland (Figures 9 and 10). High NPP regions were focused in the eastern and southern portions of the mainland, mainly in southwestern Hainan Province, Fujian Province, and most of Taiwan Province. These regions have low latitudes, sufficient heat, and abundant precipitation, which are conducive to the accumulation of organic compounds in vegetation. Low NPP regions were clustered in the western and northern parts of the continent, mainly in China's Inner Mongolia Autonomous Region, Xinjiang Uygur Autonomous Region, Gansu Province, Qinghai Province, Ningxia Hui Autonomous Region, and Tibet Autonomous Region. These areas have an arid climate, poorer soils, and low vegetation cover, which are not beneficial for the accumulation of organic material in vegetation.
The number of NPP values in the domain of 0-150 g C/m 2 was the highest, with 56.03%, followed by 300-450 g C/m 2 , 450-600 g C/m 2 , and 150-300 g C/m 2 with 22.82%, 12.90%, and 6.36%, respectively, and finally >600 g C/m 2 with the lowest percentage of 1.89%. In terms of the statistical histogram by province (Figure 9), Taiwan Province had the largest mean multi-annual NPP values of 984.44 g C/m 2 , followed by Hainan Province (724.04 g C/m 2 ) and Fujian Province (638.46 g C/m 2 ). The mean multi-year NPP value in Xinjiang Uygur Autonomous Region was the smallest, at 30.79 g C/m 2 .  The number of NPP values in the domain of 0-150 g C/m 2 was the highest, with 56.03%, followed by 300-450 g C/m 2 , 450-600 g C/m 2 , and 150-300 g C/m 2 with 22.82%, 12.90%, and 6.36%, respectively, and finally >600 g C/m 2 with the lowest percentage of 1.89%. In terms of the statistical histogram by province (Figure 9), Taiwan Province had the largest mean multi-annual NPP values of 984.44 g C/m 2 , followed by Hainan Province (724.04 g C/m 2 ) and Fujian Province (638.46 g C/m 2 ). The mean multi-year NPP value in Xinjiang Uygur Autonomous Region was the smallest, at 30.79 g C/m 2 .

Spatial Variation Characteristics
The proportion of areas with increasing and decreasing NPP trends was 72.96% and 27.04%, respectively ( Figure 11). This study presents the spatial distribution of NPP for the two subregions for 2005, 2010, 2015, and 2020, as well as the overall trend distribution for these twenty years ( Figure 11). It is obvious that the NPP in subregion A gradually increases and the NPP in subregion B progressively decreases as time passes. The significance tests of the NPP trend results indicated that the percentage of significant and very significant changes in regional NPP were 17.10% and 30.10%, respectively (Figure 12a). The results of the overlay analysis of the NPP trend and significance show a significant increase and a very significant increase of 12.07% and 28.86% in these areas, respectively (Figure 12b), while the regions with significant and very significant decreases are 3.46% and 3.31%, respectively.

Spatial Variation Characteristics
The proportion of areas with increasing and decreasing NPP trends was 72.96% and 27.04%, respectively ( Figure 11). This study presents the spatial distribution of NPP for the two subregions for 2005, 2010, 2015, and 2020, as well as the overall trend distribution for these twenty years ( Figure 11). It is obvious that the NPP in subregion A gradually increases and the NPP in subregion B progressively decreases as time passes. The significance tests of the NPP trend results indicated that the percentage of significant and very significant changes in regional NPP were 17.10% and 30.10%, respectively (Figure 12a). The results of the overlay analysis of the NPP trend and significance show a significant increase and a very significant increase of 12.07% and 28.86% in these areas, respectively (Figure 12b), while the regions with significant and very significant decreases are 3.46% and 3.31%, respectively.
There were significant and very significant increases in NPP concentrated in the Northeast Plain, Loess Plateau, and Sichuan Basin in a band-like distribution, which is related to the national policy guidance. For example, major items such as the return of cultivated land to forest and grass, the "Three-North" Shelterbelt Program, and the amelioration of saline-alkaline soil all achieved certain results [41]. Regions with significant and very significant decreases in NPP were primarily in the eastern Xinjiang Uygur Autonomous Region, western Inner Mongolia, northwestern Qinghai and Gansu Provinces, and most of the Tibet Autonomous Region, which is climate-related. These areas are located in the interior of China, with low precipitation, low vegetation cover, and vast desert areas (Taklamakan Desert, 357,300 km 2 ; Gurbantunggut Desert, 56,800 km 2 ; Badain Jaran Desert, 55,000 km 2 ; and Tengger Desert, 41,900 km 2 ), thus, the vegetation productivity level is considered low. Figure 11. Theil-Sen medium trend distribution. Figure 11. Theil-Sen medium trend distribution.   There were significant and very significant increases in NPP concentrated in the Northeast Plain, Loess Plateau, and Sichuan Basin in a band-like distribution, which is related to the national policy guidance. For example, major items such as the return of cultivated land to forest and grass, the "Three-North" Shelterbelt Program, and the amelioration of saline-alkaline soil all achieved certain results [41]. Regions with significant and very significant decreases in NPP were primarily in the eastern Xinjiang Uygur Autonomous Region, western Inner Mongolia, northwestern Qinghai and Gansu Provinces, and most of the Tibet Autonomous Region, which is climate-related. These areas are located in the interior of China, with low precipitation, low vegetation cover, and vast desert areas (Taklamakan Desert, 357,300 km 2 ; Gurbantunggut Desert, 56,800 km 2 ; Badain Jaran Desert, 55,000 km 2 ; and Tengger Desert, 41,900 km 2 ), thus, the vegetation productivity level is considered low. The impacts of temperature, precipitation, and radiation on vegetation NPP were analyzed separately, eliminating the impacts of additional elements. The main positive correlation was among vegetation NPP and temperature, precipitation, and radiation in the research region (Figure 14d). The areas where NPP showed a positive correlation with the temperature accounted for 62.19% of the total, primarily in East China and South China, the southern part of Central China, the eastern and southern parts of North and Northeast China, the northern and eastern parts of Northwest China, and the western and eastern parts of Southwest China (Figure 14a). The areas where NPP was positively correlated with precipitation were as high as 79.65%, mainly in Northeast and North China, the northern part of East China, and the northern and eastern parts of Northwest China (Figure 14b  The part of the research region in which NPP was positively correlated with precipitation, temperature, and radiation for only twenty years was larger than the area occupied by a negative correlation (Figure 15). The percentage of area with a significant positive correlation between temperature and NPP was 10.15%, and the percentage of area with a significant negative correlation was 3.53%. The precipitation and NPP had a significant positive correlation between the areas with a percentage of 31.82%, and only 1.34% of the area had a significant negative correlation. The percentage of areas with a significant positive correlation between the radiation and NPP was 13.79%, and the percentage of areas with a significant negative correlation was 1.82%.     The combined investigation of the impacts of the temperature and precipitation factors on the vegetation NPP demonstrated that the multiple correlation coefficients of temperature, precipitation, and vegetation NPP ranged from 0 to 1 (Figure 16a). The correlation between temperature, precipitation, and vegetation NPP was predominantly non-significant by the F-significance test (Figure 16b). The relative role of temperature and precipitation on NPP was further explored (Figure 17), with climatic factors affecting up to 64.18% of the area. The highest percentage of regions with a weak common influence was 42.23%. This was followed by precipitation and temperature, with 16.67% and 4.02%, respectively. The strong common influence saw the lowest percentage of regions, with only 1.25%.

Relative Contributions of Climate and Human Activities to NPP Changes
The predicted NPP change rate over the past two decades in the study area ranged from −88.5608 to 127.417, and the anthropogenic NPP slope ranged from −99.4355 to 92.722 ( Figure 18). The combination of the results of climate change and human activity contributions to NPP can be found in Figure 19. The percentage of the area where NPP increased owing to climate change was 61.87%, of which 6.63% was contributed by >80%, primarily in the eastern part of Northeast and North China. The decrease was 38.13%, mainly in the southeastern part of the Qinghai-Tibet Plateau. The region of NPP increased by human activities was 73.92%, of which 38.01% was contributed by >80%, mainly in Central China and South China, central parts of North China, north and south parts of Northeast China, east and south parts of Northwest China, and the Qinghai-Tibet Plateau. The decrease accounted for 26.08%, mainly in the northeastern part of Northwest China and the western part of North China.

Relative Contributions of Climate and Human Activities to NPP Changes
The predicted NPP change rate over the past two decades in the study area ranged from −88.5608 to 127.417, and the anthropogenic NPP slope ranged from −99.4355 to 92.722 ( Figure 18). The combination of the results of climate change and human activity contributions to NPP can be found in Figure 19. The percentage of the area where NPP increased owing to climate change was 61.87%, of which 6.63% was contributed by >80%, primarily in the eastern part of Northeast and North China. The decrease was 38.13%, mainly in the southeastern part of the Qinghai-Tibet Plateau. The region of NPP increased by human activities was 73.92%, of which 38.01% was contributed by >80%, mainly in Central China

Relative Contributions of Climate and Human Activities to NPP Changes
The predicted NPP change rate over the past two decades in the study area ranged from −88.5608 to 127.417, and the anthropogenic NPP slope ranged from −99.4355 to 92.722 ( Figure 18). The combination of the results of climate change and human activity contributions to NPP can be found in Figure 19. The percentage of the area where NPP increased owing to climate change was 61.87%, of which 6.63% was contributed by >80%, primarily in the eastern part of Northeast and North China. The decrease was 38.

Relationship between Drivers and NPP Changes at Each Stage
The driving factors affecting vegetation NPP changes may differ at different stages, and the main contributors may differ for each time period. In this study, the relationships between each driver and vegetation NPP changes were analyzed using a structural equation PLS-SEM with a five-year cycle by integrating soil, climate, socioeconomic, and topographical factors.
This research first validated the reliability and validity of the PLS-SEM model ( Table  5). VIF was used to detect multicollinearity in the variables and provided values between 1 and 5, indicating that there was no significant covariance between the factors. A composite reliability (CR) greater than 0.7 and not lower than 0.6 indicated that the model was reliable, and an average variance extracted (AVE) greater than 0.5 demonstrated that the model was valid. The data in this study satisfied the above conditions well and prove that the model was reasonable and reliable [42].

Relationship between Drivers and NPP Changes at Each Stage
The driving factors affecting vegetation NPP changes may differ at different stages, and the main contributors may differ for each time period. In this study, the relationships between each driver and vegetation NPP changes were analyzed using a structural equation PLS-SEM with a five-year cycle by integrating soil, climate, socioeconomic, and topographical factors.
This research first validated the reliability and validity of the PLS-SEM model ( Table  5). VIF was used to detect multicollinearity in the variables and provided values between 1 and 5, indicating that there was no significant covariance between the factors. A composite reliability (CR) greater than 0.7 and not lower than 0.6 indicated that the model was reliable, and an average variance extracted (AVE) greater than 0.5 demonstrated that the model was valid. The data in this study satisfied the above conditions well and prove that the model was reasonable and reliable [42].

Relationship between Drivers and NPP Changes at Each Stage
The driving factors affecting vegetation NPP changes may differ at different stages, and the main contributors may differ for each time period. In this study, the relationships between each driver and vegetation NPP changes were analyzed using a structural equation PLS-SEM with a five-year cycle by integrating soil, climate, socioeconomic, and topographical factors.
This research first validated the reliability and validity of the PLS-SEM model (Table 5). VIF was used to detect multicollinearity in the variables and provided values between 1 and 5, indicating that there was no significant covariance between the factors. A composite reliability (CR) greater than 0.7 and not lower than 0.6 indicated that the model was reliable, and an average variance extracted (AVE) greater than 0.5 demonstrated that the model was valid. The data in this study satisfied the above conditions well and prove that the model was reasonable and reliable [42].
The PLS-SEM model of the relationship between the latent and explicit variables for the changes in NPP in 2001-2020 is shown in Figure 20

Comparison with Other NPP Reports
The data sources used in NPP studies can be broadly divided into two categories: existing NPP products [43][44][45][46][47] (Table 6), which have a fixed spatial resolution, and estimations of NPP using NDVI and FPAR. To analyze the relationship between NPP changes and influencing factors, studies have mainly explored the overall response mechanism of NPP to climate change and human activities; in addition, they have considered mainly the effects of temperature and precipitation [36,43,46,48,49]. Other studies have considered additional influencing factors, such as relative humidity, sunshine hours, chlorophyll-a, photosynthetic availability, wind stress curl, salinity, nitrates, phosphates, and CO2 concentration [44,47,50]. Compared with other studies, this paper conducted a relatively comprehensive study of the climatic and anthropogenic influences on vegetation NPP; however, other influencing factors, such as CO2 concentrations, nitrates, and phosphates, could be considered in subsequent studies. The amount of fertilizer used has an effect on vegetation NPP changes; unfortunately, raster data on fertilizer inputs such as N and P were lacking in this study. In a forthcoming follow-up study, we intend to obtain raster data on various types of fertilizer inputs to enhance the accuracy and reliability of the influence factors of vegetation NPP.

Comparison with Other NPP Reports
The data sources used in NPP studies can be broadly divided into two categories: existing NPP products [43][44][45][46][47] (Table 6), which have a fixed spatial resolution, and estimations of NPP using NDVI and FPAR. To analyze the relationship between NPP changes and influencing factors, studies have mainly explored the overall response mechanism of NPP to climate change and human activities; in addition, they have considered mainly the effects of temperature and precipitation [36,43,46,48,49]. Other studies have considered additional influencing factors, such as relative humidity, sunshine hours, chlorophyll-a, photosynthetic availability, wind stress curl, salinity, nitrates, phosphates, and CO 2 concentration [44,47,50]. Compared with other studies, this paper conducted a relatively comprehensive study of the climatic and anthropogenic influences on vegetation NPP; however, other influencing factors, such as CO 2 concentrations, nitrates, and phosphates, could be considered in subsequent studies.

Limitations and Prospects
Although the accuracy of kNDVI estimates at the site scale was higher than that of NDVI, unfortunately, ChinaFLUX shared flux observations only for certain stations. This makes it difficult to rely on observational data alone when assessing the quality of NPP data at the national scale. In the future, more measured data will be sought as much as possible to improve the quality of the estimated NPP data.
For accuracy validation, direct real NPP data were not available in this research, and were obtained by converting flux data. It is worth noting that the value of the NPP/GPP ratio α used in conversion varies with geography, meteorology, and vegetation type, while in this study we used the regional average NPP/GPP ratio (α = 0.64), which could potentially have an impact on the NPP conversion value.

Conclusions
Based on an improved CASA model driven by the new kNDVI vegetation index, this paper estimated the NPP of China's terrestrial ecosystem vegetation from 2001 to 2020 and used the Theil-Sen median, Mann-Kendall trend, second-order partial correlation, multiple correlation analysis, differential equation method, and PLS-SEM method to analyze the temporal and spatial evolution of vegetation NPP in China's terrestrial ecosystems. The results show that kNDVI can effectively improve the saturation problem, cope with complex phenological changes, and improve the precision of NPP estimation. The NPP of Chinese vegetation presented a fluctuating upward trend on the time scale from 2001 to 2020. It increased by a mean of 0.7235 g C/m 2 each year, showing a low northwest to high southeast pattern in the spatial distribution, with the highest average multi-year NPP value by province in Taiwan Province and the lowest in Xinjiang Uygur Autonomous Region. The changes in NPP were dominated by an increasing trend (72.96%), and the areas showing significant and very significant increases were concentrated in the Northeast Plain, Loess Plateau, and Sichuan Basin. In Northeast and North China, vegetation NPP was mainly influenced by precipitation, while in Southwest and South China it was mainly influenced by temperature and in Central China by solar radiation. With the passage of time, the dominant factors influencing the changes in NPP change and have different impacts. The main influencing factors of the climate in 2001-2005 were evaporation and precipitation, and the main influencing factors of the climate in 2005-2010 were evaporation and temperature; moreover, the influences of both phases were positive. However, in 2010-2015 the impact of the climate was negative. Climate change and human activities have a major positive effect on NPP change. The areas where the NPP increased due to climate change were mainly located in Northeast and North China, and the areas where the NPP increased due to human activities were mainly located in Central and South China. The results of this research can be used in China and around the world to better prepare for and respond to climate change.