A 33-Year NPP Monitoring Study in Southwest China by the Fusion of Multi-Source Remote Sensing and Station Data

Abstract

Knowledge of regional net primary productivity (NPP) is important for the systematic understanding of the global carbon cycle. In this study, multi-source data were employed to conduct a regional NPP study in southwest China, with a 33-year time span and a 1-km scale. A multi-sensor fusion framework was applied to obtain a new normalized difference vegetation index (NDVI) time series from 1982 to 2014, combining the advantages of different remote sensing datasets. As another key parameter for NPP modeling, the total solar radiation was calculated utilizing the improved Yang hybrid model (YHM), based on meteorological station data. The accuracy of the data processes is proved reliable by verification experiments. Moreover, NPP estimated by fused NDVI shows an obvious improved accuracy than that based on the original data. The spatio-temporal analysis results indicated that 67% of the study area showed an increasing NPP trend over the past three decades. The correlation between NPP and precipitation was significant heterogeneous at the monthly scale; specifically, the correlation is negative in the growing season and positive in the dry season. Meanwhile, the lagged positive correlation in the growing season and no lag in the dry season indicated the important impacts of precipitation on NPP. What is more, we found that there are three distinct stages during the variation of NPP, which were driven by different climatic factors. Significant climate warming led to a great increase of NPP from 1992 to 2002, while NPP clearly decreased during 1982–1992 and 2002–2014 due to the frequent droughts caused by the precipitation decrease.

1. Introduction

As a key component of the global carbon cycle, the terrestrial ecosystem is the main force that can uptake free carbon from the atmosphere and convert it into organic compounds [1,2]. An improved understanding of the terrestrial ecosystem carbon cycle is important under the background of intense global climate change [3,4,5]. Net primary productivity (NPP), which is the residual amount of organic matter produced by vegetation photosynthesis minus its autotrophic respiration consumption, is an important ecological indicator for the status of a terrestrial ecosystem carbon budget [6,7]. NPP can be precisely acquired by field measurements at a site level, but it is not feasible for regional or larger scales, as it costs a lot and is unable to obtain spatially full coverage and continuously long time series. In this case, model-based estimation is an efficient approach, and a large number of models have been proposed in previous studies [8]. With the development of remote sensing, satellite data-based models have been extensively applied to terrestrial NPP researches. The reason for the popularity of the satellite data-based models is that remote sensing images can provide continuous, dynamic, and comprehensive land-surface information for any region around the earth [9,10,11]. Satellite land-cover data and spectral vegetation index products (i.e., the normalized difference vegetation index, NDVI) are the most commonly used core data when modeling NPP of a large region [8].

Numerous studies have explored the spatio-temporal patterns of NPP based on moderate-resolution satellite NDVI datasets [12,13,14], such as the Advanced Very High Resolution Radiometer (AVHRR) data from the National Oceanic and Atmospheric Administration (NOAA) satellites and the Moderate Resolution Imaging Spectroradiometer (MODIS) data from the Terra/Aqua satellites [15,16]. However, the related studies have mostly concentrated on the variation of NPP at global or continental scales [13,17,18,19]. Although some research has concerned the regional carbon cycle [20,21,22], the in-depth studies are still lacking. Nevertheless, regional terrestrial ecosystem carbon budgets are important, because the characteristics of NPP can highly vary in space and time due to different natural environments and human landscapes [23,24]. Further study of regional carbon cycles is therefore necessary for the systematic understanding of the global carbon cycle, since an integrated and dense carbon observation and analysis system is urgently required [25].

In regional carbon cycle studies with satellite data-based models, remote sensing data quantity and quality are the determinants for a significant and in-depth analysis. MODIS datasets have been applied in many studies of the spatio-temporal variations of regional NPP and its potential causal factors, at spatial resolutions of 250 m to 1 km [24,26,27,28]. Although these studies have captured the spatial information of NPP, they have been unable to analyze the NPP variation and its relationship with environmental factors in the long run, because only data covering the last decade are available, and no data are available before the year 2000. The same limitation also exists in the studies with Satellite Pour l’Observation de la Terre (SPOT) Vegetation products [29,30,31]. Studies using Landsat data can offer even more detailed spatial information, but it is difficult to acquire a continuous long-term series, because of its poor temporal resolution and the influence of cloud cover [32,33,34]. Since the 1980s, AVHRR datasets have been extensively employed to study the regional NPP of a long time period [20,35,36,37], but they have been generally used for much larger regions, due to the coarse spatial resolution of 8 km. Furthermore, it has been proven in many papers that a coarse resolution can lead to an obvious accuracy loss when modeling NPP, as a result of the spatial heterogeneity of the data [38,39,40]. In general, further study of the regional carbon cycle is limited by the inter-inhibitive characteristics of the different sensors. In this condition, integrating remote sensing data from different sensors would be an efficient way to settle the problem [41,42,43]. Related attempts have been made to relieve the trade-off between the temporal and spatial resolutions in NPP study [44,45,46], which proved the feasibility of applying multi-source remote sensing data. However, the insufficient characteristics of spatial resolution and time span have not been solved previously, which would be more important for NPP studies. Thus, it is of great significance to obtain a long-term NPP series with suitable resolution by combing multi-source data.

Yunnan province, located in the most southwestern part of China, is one of the most important carbon sinks in the continent [47]. However, this region has suffered from more and more frequent droughts in recent decades [48,49]. A continuous four-year extreme drought occurred from 2009 to 2012, which severely impacted up to 56% of the region [50]. Although a number of studies have investigated the drought impacts in recent years [50,51], meritorious research into the long-term variation of the regional carbon cycle and its relationship with climatic factors is extremely rare. As a result, to address these issues, the main objectives of this paper are: (1) to develop an innovative framework for the generation of a continuous 33-year NPP time series for Yunnan province at a 1-km spatial resolution with the Carnegie-Ames-Stanford Approach (CASA) model by fusing multi-source remote sensing data and station data; and (2) to carry out an in-depth analysis of the spatial and temporal characteristics of the terrestrial ecosystem carbon cycle in the study area, as well as its relationships with the climate.

2. Data Sources

2.1. Study Area

Yunnan province, which is located in southwest China between 21.13°–29.25°N and 97.52°–106.18°E, was chosen as the study area in this research (Figure 1). Yunnan covers a total area of 394,000 km², and mountainous landforms occupy more than 90% of the region. The unique mountainous landforms lead to its diverse climate characteristic and provide conditions for the growth of many different vegetation species. Vegetation covers most of the area (approximately 94%), and the rich and varied vegetation types range from tropical species to frigid species [52]. Most of the territory is located in the subtropical or tropical zones, with the Tropic of Cancer running through its southern part. The synergistic effects of the tropical/subtropical climate, the monsoon climate, and the mountain climate lead to highly complex climate patterns in the area. Although the temperature is moderate and rainfall is abundant, the uneven intra-annual allocation usually results in wet summers and dry winters [50].

2.2. Data Sources

2.2.1. NDVI Datasets

The Global Inventory Modelling and Mapping Studies 3rd generation (GIMMS3g) NDVI product from the NOAA/AVHRR and MODIS monthly NDVI collection (MOD13A3) were selected as the basis to composite the 33-year NDVI time series at a 1-km scale. The GIMMS3g dataset from 1982 to 2012 was obtained from the National Aeronautics and Space Administration (NASA) Ames Ecological Forecasting Lab. The spatial resolution of this dataset is 8 km, and the time interval is half a month. The dataset has been proven to have a better availability and quality than other AVHRR-based NDVI products [53]. The MOD13A3 collection from 2000 to 2014 (no data are available for January 2000) was acquired from the NASA Earth Observing System (EOS) program, at a spatial resolution of 1 km. The data, which are obtained based on the spectral bands (red band and near-infrared band) that are primarily designed for the study of vegetation and the land surface [54], have been widely applied in numerous vegetation studies [14]. In this study, the maximum value compositing (MVC) technique was employed to obtain the monthly GIMMS3g NDVI, to match the time interval of the MODIS data.

2.2.2. Meteorological Datasets

The meteorological datasets over the study period of 1982–2014 were obtained from the China Meteorological Administration (CMA), including monthly and daily precipitation and air temperature, and daily surface pressure, air relative humidity, and sunshine duration from the 29 uniformly distributed meteorological stations (Figure 1b). In addition, radiation datasets of the five specified radiation stations (Figure 1b) were also compiled. The station records were carefully interpolated into the same spatial resolution as the MODIS NDVI (1 × 1 km), using the Australian National University SPLINe (ANUSPLIN) package [55], with elevation or slope data as the independent covariates. This package can be used to undertake professional meteorological interpolation using the thin plate smoothing splines surface-fitting technique. The use of independent covariates can further improve the precision. After repeated attempts, it was concluded that elevation was the optimal covariate for temperature, and slope was the optimal covariate for precipitation and radiation.

2.2.3. Other Data

The WESTDC2.0 land-cover map was derived by the Chinese Academy of Sciences (CAS) Environmental and Ecological Science Data Center for West China (WESTDC) [56]. The map takes full advantage of the 1:100,000 land resources data surveyed by CAS, and is integrated with multi-source satellite classification information. In this study, the data were synthesized into eight classes: evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), needle-leaf forest (NF), mixed forest (MF), shrub, grass, crop, and other land covers.

The measurement-based biomass/NPP datasets from Luo’s study [57], which have been used in many studies [51,58], were employed as the validation data in this study. All the records over the study area were from the Yunnan Ministry of Forestry for the year of 1983. The data include the forest biomass/NPP for most of the plant components, and the location and dominant species of each site. As the records of NPP were provided with the unit of dry matter (t DM ha⁻¹ year⁻¹), a conversion factor of 50 was needed to change this into carbon content (gC m⁻² year⁻¹) [58,59].

3. Method

3.1. CASA Model

The Carnegie-Ames-Stanford Approach (CASA) model, which was developed on the basis of light-use efficiency, served to estimate the monthly NPP in the study area [60,61]. The calculation of NPP can be expressed as the product of absorbed photosynthetic active radiation ($APAR$, MJ m⁻²) and the light-use efficiency ($\epsilon$, gC MJ⁻¹), as follows:

$$NPP\left(x,t\right)=APAR\left(x,t\right)\times \epsilon \left(x,t\right),$$

(1)

where $NPP\left(x,t\right)$ is the total fixed NPP of pixel $x$ in month $t$, $APAR\left(x,t\right)$ is the total amount of absorbed photosynthetic active radiation over the period, and $\epsilon \left(x,t\right)$ is the actual light-use efficiency. The calculation of $APAR$ and $\epsilon$ is shown below:

$$APAR\left(x,t\right)=R_{\mathrm{s}}\left(x,t\right)\times 0.5\times FPAR\left(x,t\right),$$

(2)

$$\epsilon \left(x,t\right)={\epsilon }^{\ast }(x,t)\times T_1(x,t)\times T_2(x,t)\times W(x,t),$$

(3)

where $R_{\mathrm{s}}\left(x,t\right)$ is the total solar radiation of pixel $x$ in month $t$; the coefficient 0.5 is the approximate ratio of photosynthetic active radiation (0.4–0.7 μm) to total solar radiation; and $FPAR\left(x,t\right)$ is the fraction of photosynthetic active radiation absorbed by the vegetation canopy, which is determined by the NDVI and vegetation types:

$$FPAR(x,t)=\min [\left(SR(x,t)-S{R}_{\min }\right)/\left(S{R}_{\max }-S{R}_{\min }\right),0.95],$$

(4)

$$SR(x,t)=[1+NDVI(x,t)]/[1-NDVI(x,t)],$$

(5)

where $SR\left(x,t\right)$ is the simple ratio of NDVI; $S{R}_{\max }$ and $S{R}_{\min }$ are the constants related to the vegetation type, which is shown in

Table 1. The value of maximum light utilization efficiency (${\epsilon }^{\ast }$), $S{R}_{\min }$, $S{R}_{\max }$ for different vegetation types.

Vegetation Type	EBF	DBF	NF	MF	Shrub	Grass	Crop
${\epsilon }^{\ast }$ (gC MJ−1)	0.985	0.692	0.485	0.768	0.429	0.542	0.542
$S{R}_{\min }$	1.050	1.050	1.050	1.050	1.050	1.050	1.050
$S{R}_{\max }$	5.170	6.910	6.630	4.670	4.490	4.460	4.460

Note: ${\epsilon }^{\ast }$ is the maximum light-use efficiency; $S{R}_{\min }$ is the factor $SR$ for unvegetated land areas; $S{R}_{\max }$ approximates the values of $SR$ when all solar radiation is intercepted; EBF: evergreen broadleaf forest; DBF: deciduous broadleaf forest; NF: needle-leaf forest; MF: mixed forest.

[62].

In Equation (3), ${\epsilon }^{\ast }\left(x,t\right)$ is the maximum light-use efficiency, the value of which varies with the vegetation type according to the previous study of ecosystems in China (Table 1) [62]; $T_1\left(x,t\right)$ and $T_2\left(x,t\right)$ are the temperature stress factors; and $W\left(x,t\right)$ is the moisture stress factor:

$$T_1(x,t)=0.8+0.02\times T_{\mathrm{opt}}(x)-0.0005\times T_{\mathrm{opt}}(x)\times T_{\mathrm{opt}}(x),$$

(6)

$$T_2(x,t)=1.1814/\{1+e^{[0.2(T_{\mathrm{opt}}(x)-10-T(x,t))]}\}/\{1+e^{[0.3(-T_{\mathrm{opt}}(x)-10+T(x,t))]}\},$$

(7)

$$W(x,t)=0.5+0.5EET\left(x,t\right)/PET\left(x,t\right),$$

(8)

where $T_{\mathrm{opt}}\left(x\right)$ is the temperature when NDVI reaches its maximum for the year; $T\left(x,t\right)$ is the monthly mean temperature; $EET\left(x,t\right)$ and $PET\left(x,t\right)$ are the soil properties derived from the sub-model of regional evapotranspiration [63].

The spatial and temporal resolution of the NPP from the CASA model is determined by the resolution of the NDVI. Thus, in Section 3.2, image processing algorithms are introduced to integrate the respective advantages of the MODIS and GIMMS3g datasets, to obtain a long-term NDVI time series with a suitable resolution. In addition, in Section 3.3, the total solar radiation is precisely calculated by the improved YHM model. The overall workflow of NPP estimation is depicted in Figure 2.

3.2. Generation of the 33-Year NDVI Time Series at a 1-km Scale

The three processes described below were adopted to eliminate the problems with regard to data quality, sensor differences, and the coarse spatial resolution, for the NDVI datasets.

3.2.1. NDVI Filtering

The application of NDVI time-series data is usually limited by the existence of unwanted noise and errors caused by the cloud presence and other atmospheric contamination [64]. In order to obtain a high-quality NDVI time series, the moving weighted harmonic analysis (MWHA) [65] method was employed to correct these contaminated values in the GIMMS3g and MODIS NDVI datasets. The algorithm was proposed based on a modification of the harmonic analysis method, which has been proved to be a better strategy for different NDVI datasets with various time intervals, and is reliable for pixels in the vegetation dormancy period [65]. The effects of the filter process are shown in Figure 3, where it can be observed that the sudden abnormal drops in the NDVI time series have been exactly adjusted. After the filtering, the time series show much more similar variation curves approaching the actual vegetation variation.

3.2.2. Normalization

A large gap between the NDVI of the two sensors for the specific pixel can be observed in Figure 3, and the NDVI values of MODIS are generally lower than GIMMS3g, with a mean difference of 0.07 in the study area. This is likely to be caused by the disparity in the characteristics of the two sensors, and the fact that no atmospheric correction has been applied to the GIMMS3g data [66,67]. Fensholt and Proud concluded that the temporal trends derived from the GIMMS NDVI agree well with the MODIS data [68], overall, so a unary linear regression normalization model can be used to express the relationship between the two datasets. In this study, the GIMMS3g data were normalized to be consistent with MODIS using a pixel-by-pixel linear regression method, at the original spatial resolution of the GIMMS3g data (8 km). For each pixel, a linear relationship was obtained based on the data in the mutual time period from February 2000 to December 2012 (155 pairwise images). Then it was applied to the GIMMS3g data for the corresponding pixels before the year 2000. Thus, the data time series from 1982 to 1999 can be considered as MODIS-like NDVI, with data values consistent with MODIS, but with the same resolution as GIMMS3g.

3.2.3. Multi-Sensor Fusion

Although the sensor differences were removed, gaps still existed in the spatial resolution between the MODIS NDVI and the obtained MODIS-like time series. On account of the 8-km resolution being too coarse, which can lead to overestimation when modeling NPP [38,39], multi-sensor fusion was an effective way to improve the spatial resolution of the MODIS-like data [69,70,71]. Many vegetation studies have applied the spatial and temporal adaptive reflectance fusion model (STARFM) and the extended STARFM (ESTARFM) to vegetation index fusion and prediction [41,72,73,74,75]. In this study, a spatio-temporal information fusion method based on a non-local means filter was employed to improve the spatial resolution of the MODIS-like data [76].

Multi-sensor fusion can be used to predict the MODIS NDVI value at $t_1$ based on the MODIS-like data at $t_1$ and the reference MODIS and MODIS-like images acquired at $t_0$. The prediction of the fine-resolution NDVI before the year 2000 can be expressed as:

$$F\left(x_{p/2},y_{p/2},t_1\right)={\displaystyle \sum _{i=1}^p{\omega }_i\times \left(F\left(x_i,y_i,t_0\right)+C\left(x_i,y_i,t_1\right)-C\left(x_i,y_i,t_0\right)\right)},$$

(9)

where $F$ and $C$ represent the fine-resolution and coarse-resolution NDVI, respectively; $t_0$ is the acquisition date of the reference data; $t_1$ is the prediction date; $\left(x_{p/2},y_{p/2}\right)$ is the location of the predicted pixel; $\left(x_i,y_i\right)$ denote the pixel location; $p$ is the size of the moving window; and ${\omega }_i$ is the spatial weighting function. The innovation of the applied fusion algorithm is the more reasonable calculation of ${\omega }_i$, which takes full consideration of the spatial relationship between pixels based on the concept of the non-local means filter. Considering the seasonal inconsistency of NDVI patterns, the fine-resolution NDVI of each month from 1982 to 1999 was predicted referring to the MODIS data for the corresponding month in the nearest year.

3.3. Accurate Calculation of Total Solar Radiation with the Improved YHM Model

The surface total solar radiation ($R_{\mathrm{s}}$) is an indispensable parameter of the CASA model. However, the sparse distribution of the $R_{\mathrm{s}}$ observation stations (five in the study area) resulted in it being difficult to interpolate accurate raster radiation data [77]. Therefore, the improved Yang hybrid model (YHM) was employed to calculate the $R_{\mathrm{s}}$ of the 29 climatological stations in Yunnan province [77], using the daily climatological records. The model considers the attenuation from each atmospheric component when solar radiation passes through the atmosphere. The model can be described as follows:

$$R_{\mathrm{s}}={\tau }_c{\displaystyle \underset{\Delta t}{\int }({\tau }_{b,clear}+{\tau }_{d,clear})}{I}_0\mathrm{d}t,$$

(10)

where $I_0$ is the solar irradiance at the top of atmosphere; $\Delta t$ is the time period of the calculated solar radiation; ${\tau }_{b,clear}$ and ${\tau }_{d,clear}$ are the solar beam radiative transmittance and the solar diffuse radiative transmittance under clear skies, respectively, which were determined by the transmittances of ozone, water vapor, gas mixture, the Rayleigh effect, and aerosols in the atmosphere; and ${\tau }_c$ is the radiative transmittance of cloud, which is a linear function of sunshine duration. The improved estimation of parameter ${\tau }_c$ was undertaken according to the method proposed by Wang et al. [78].

4. Results

In this part, the reliability of estimated NPP was firstly verified in Section 4.1, as well as the accuracy of NDVI fusion and total solar radiation calculation. Then, in Section 4.2, the spatial and temporal variation of NPP in Yunnan was characterized. In order to unveil the climatic impacts on NPP, their correlations were analyzed at an annual and a monthly scale, respectively in Section 4.3.1 and Section 4.3.2 Then the lagged effects of precipitation were also considered in Section 4.3.3. Finally, in the Section 4.4, the dominating climate factor for NPP variation was found for different time stages in the past three decades.

4.1. Results Validation

4.1.1. Simulated Validation of the Multi-Sensor Fusion

In order to validate the feasibility of the multi-sensor fusion process in long-term NDVI prediction, simulated experiments were conducted based on the MODIS and GIMMS3g datasets for the mutual time period of 2000 to 2012. For each month, the fine-resolution NDVI ($NDV{I}_{\mathrm{fusion}}$) from 2000 to 2011 was simulated with reference to the MODIS and MODIS-like images for the corresponding month in the year 2012. The predicted fusion data ($NDV{I}_{\mathrm{fusion}}$) were then validated with the true MODIS NDVI, both qualitatively and quantitatively. The statistics of root-mean-square error (RMSE) and mean absolute difference (MAD) for each year were computed and are shown in Figure 4. The fused NDVI presents an absolute error of less than 0.06 compared with the original MODIS data, and the RMSE is around 0.07. Most importantly, the fine fusion results maintain a stable accuracy as the years between the prediction date and reference date increase. It is, therefore, reasonable to predict the fine-resolution monthly NDVI before the year of 2000 using the multi-sensor fusion method. What is more, a qualitative comparison of the results for April 2000 is shown in Figure 5 as an example, where it can be observed that the spatial resolution of the predicted NDVI has been clearly improved. The spatial distribution of the data is highly consistent, and only slight differences can be found. As shown in the two green rectangles in Figure 5e,f, although some regions with low NDVI have been overestimated in the fusion result, the vegetated land area is well predicted. It means that the fusion results are useable in the vegetation related studies.

4.1.2. Cross-Validation of the Total Solar Radiation Calculation

Cross-validation experiments were applied to examine the accuracy of the YHM model based on the quality-controlled monthly solar radiation measurements ($R_{\mathrm{station}}$) [79], of five stations in Yunnan (station numbers 56651, 56739, 56778, 56959 and 56985, respectively). The records of four stations were used to obtain parameter ${\tau }_c$ and the model solar radiation ($R_{\mathrm{model}}$), and then the $R_{\mathrm{station}}$ records of the remaining one station (test station) were applied to validate the calculated $R_{\mathrm{model}}$ of the corresponding site. The experiment was undertaken five times until every station was tested. The statistics of the correlation coefficient (r), RMSE, and mean absolute relative difference (MARD, the mean absolute value of the ratio between the error and true data) are listed in

Table 2. Statistics for the cross-validation of the improved YHM model.

Test Station No.	r	RMSE (MJ m−2)	MARD (%)
56651	0.88	48.24	7.31
56739	0.84	52.76	8.71
56778	0.94	48.60	7.67
56959	0.85	46.86	7.96
56985	0.83	53.72	8.41

Note: Test station No. is the station number that used to validate the solar radiation calculated by the other four stations; RMSE is the root mean square error; and MARD denotes the mean absolute relative difference, which is the mean absolute value of the ratio between the error and true data.

. The results show that the modeled solar radiation has a good consistency with the observed records, with r above 0.83 and MARD at around 8%.

4.1.3. Validation of the Estimated NPP with in situ Data

The reliability of the estimated NPP products was verified using the in situ data measured by the Yunnan Ministry of Forestry [57]. Totally 59 plots were picked for validation, through matching the vegetation type with the land-cover map used in this study. As displayed in Figure 6, the estimated NPP shows fine consistency with the in situ data, with r reaching 0.79 (p < 0.001). Moreover, the good consistency also could be observed among different vegetation types, as shown in

Table 3. Statistics of the estimated NPP performances at 59 field sites.

Vegetation Types	Measured Value (gC m−2 year−1)			In This Study (gC m−2 year−1)
	Mean	Range	STD	Mean	Range	STD	MAD	Bias
EBF	996	524–1233	161	1080	675–1474	213	63	14
DBF	778	693–864	121	774	685–862	125	5	−5
NF	536	313–854	152	611	428–779	80	105	74
ALL	563	313–1233	184	633	428–1474	122	99	67

Note: STD is the standard deviation; MAD is the mean value of the absolute difference for the modelled NPP and observed data; Bias is the mean value of the modeled data minus observed data.

. The estimated NPP demonstrates the similar statistics with the in situ data, including the mean value, data range and standard deviation. The overall MAD is 99 gC m⁻² month⁻¹ and the bias is 67 gC m⁻² month⁻¹. Although there are some disagreements, the scale effect and representative errors might be the major reason for it. In general, the estimated NPP shows an overall good consistency with the in situ data, with high r and low difference, which demonstrates its usability in analyzing the NPP patterns in Yunnan province.

4.2. NPP Spatial Distribution and Variation Trends

4.2.1. NPP Spatial Distribution

The spatial distribution of the mean annual NPP for the past 33 years is shown in Figure 7a. Generally, the NPP in Yunnan province gradually decreases from the southwest to the northeast. The mean annual NPP is generally higher than 1000 gC m⁻² year⁻¹ in most of the southwest area, which is located close to the frontier, where the latitude and altitude are relatively low. The lack of human activities and the warm climate in this area benefit the growth of vegetation. However, the mean annual NPP is less than 500 gC m⁻² year⁻¹ in extensive regions of northwest Yunnan, which is a part of the Qinghai–Tibet Plateau, with elevations mostly higher than 4000 m. The cold and harsh climate in this area limits the vegetation growth, which is also the reason for the low NPP in the northeast area. There are also many regions with an NPP of less than 500 gC m⁻² year⁻¹ in central Yunnan, such as the provincial capital of the city of Kunming, as a result of urban construction and expansion. In general, most of Yunnan shows a relatively high NPP exceeding 800 gC m⁻² year⁻¹, but some regions have lower NPP values due to the harsh climate or human activities.

4.2.2. Annual NPP Variation

The inter-annual changes of the mean annual NPP in Yunnan province are shown in Figure 8a. Overall, the mean annual NPP showed fluctuating growth from 1982 to 2014, with a total increasing trend of 0.98 gC m⁻² per year (r = 0.38, p < 0.05). However, the inter-annual variation of NPP was not consistent over the entire study period, but three distinct stages were found after filtering the variation curve, with break points at 1992 and 2002. It can be clearly observed in the figure that the NPP experienced a decreasing trend from 1982 to 1992 (slope = −3.04 gC m⁻² year⁻¹); it then sharply increased at a rate of 5.70 gC m⁻² year⁻¹ until 2002; and finally slightly decreased again between 2002 and 2014 (slope = −2.22 gC m⁻² year⁻¹). Therefore, the increasing trend of the NPP over the study period was mainly due to the increment from 1992 to 2002, because the NPP in the other two stages presented completely opposite variation trends.

The spatial pattern of the annual NPP trends (the slope of the NPP inter-annual variation) in the past 33 years is shown in Figure 7b. The trends are calculated by the robust regression of iteratively reweighted least squares [80], which could eliminate the impacts of outliers. The annual NPP increased in 67.16% of the study area, which was more than twice the region with decreased NPP. In particular, 37.32% of the area showed a significantly increased NPP trend (r > 0, p < 0.05), while the rate for the significantly decreased area was only 9.48% (r < 0, p < 0.05). Almost all of southwest Yunnan showed increasing trends, with values higher than 2 gC m⁻² per year. This indicates that the NPP in southwest Yunnan is not only high, but has also presented a significant increasing trend over the past 33 years. Obvious increasing trends also occurred in large regions in the northeast, which shows a relatively low NPP. Meanwhile, there are also regions showed decreasing trends, mostly in the high-altitude district in the northwest and the urban area in central Yunnan. The hostile climate and urban expansion might be the reason for the NPP decline. It is found that 91.02% of these areas have suffered from a decreased precipitation or temperature in last three decades, especially area with decreased temperature occupied 75.75%.

4.3. Relationship between NPP and Climate

4.3.1. Correlations between NPP and Climatic Factors at an Annual Scale

The interaction between NPP and climatic factors is one of the most important issues for carbon cycle research. The inter-annual changes of mean temperature and accumulated precipitation for the study area are shown in Figure 8b,c, and are divided into the three stages with the same break points in NPP annual variation. It can be seen that the temperature steadily increased in all stages, with an overall increasing trend of 0.039 °C per year (r = 0.81, p < 0.01), and the total increment reached 1.28 °C from 1982 to 2014. The benefit of continuous climate warming was responsible for the overall growth in NPP [12]. Temperature also showed a significant positive correlation with NPP at an annual scale, with r = 0.41 (p < 0.05). Precipitation presented three completely parallel stages to NPP, with trends of −14.61 mm year⁻¹, 15.28 mm year⁻¹ and −12.56 mm year⁻¹, respectively. The relative magnitudes of the variation trends for the three stages also agreed with the NPP. However, the precipitation showed an overall variation trend that was the opposite to NPP, with a decreasing trend of −4.79 mm per year (r = −0.43, p < 0.05). An abnormally negative correlation was observed between annual NPP and precipitation, in spite of the fact that the relationship was not significant (r = −0.22, p > 0.1). This indicates that increased precipitation led to reduced NPP for vegetation in Yunnan province.

4.3.2. Correlations between NPP and Climatic Factors at a Monthly Scale

In order to conduct further studies of the relationships between NPP and climatic factors, the intra-annual variations of monthly NPP, climatic factors, and their correlations are shown in Figure 9, as well as the important parameter of total solar radiation. It can be seen that there was an obviously higher NPP in summer and low accumulation in winter, with the NPP difference as large as 60 gC m⁻² per month. This characteristic of NPP is the result of the similar uneven distribution of precipitation and temperature at a monthly scale. The mean temperature in summer was more than the double that in winter, with a difference of more than 10.46 °C. The heterogeneous intra-annual distribution of precipitation was even greater. The half-year from May to October featured 82.88% of the annual precipitation. The suitable temperature and abundant precipitation provide good conditions for vegetation growth from May to October, while the cold and dry climate suppresses photosynthesis and leads to low NPP in the other months. According to the previous vegetation study in Yunnan [81], the period from May to October is the “growing season”, with high NPP and a superior hydrothermal environment, and the other six months can be called the “dry season”, with low NPP and a harsh climate.

Figure 9 also shows the difference in the correlation coefficients (r) between monthly NPP and the climatic variables for each month. Temperature presents a relatively stable positive correlation with NPP throughout the year, with a relatively weaker impact in summer when the warm climate is at the optimum level for photosynthesis. Nevertheless, highly fluctuating correlation coefficients between NPP and precipitation and total solar radiation are observed in the intra-annual variation, and completely antipodal distributions are found between them. Total solar radiation fully controls the change of NPP from June to October, with amazingly high correlation coefficients of greater than 0.97 (p < 0.01). The relationship in May is also significant (r = 0.69, p < 0.01). In contrast, precipitation shows a significant negative impact on the vegetation growth in the six months of the growing season. This phenomenon is the reason for the negative correlation between precipitation and NPP at an annual scale. Due to NPP in the growing season amounts to 65.7% of the annual amount, the significant negative impacts of precipitation during this period might dominated the weak negative correlation at an annual scale.

The negative impact of precipitation on NPP in the growing season is the result of the warm temperature and abundant precipitation during this period, whose mean values are 18.61 °C and 200.31 mm per month. The least rainfall in the growing season occurred in the year 2011, but the average amount was still as high as 67.09 mm per month. In this condition, total solar radiation, as the energy source of photosynthesis, certainly dominates the growth of vegetation, with a very close correlation to NPP. However, due to the rainfall already being abundant for photosynthesis during this period, more precipitation cannot benefit the vegetation and actually suppresses its growth. This is because the increased precipitation means more overcast skies, which prevents the transfer of sunlight and lessens the total amount of solar radiation arriving at the vegetation canopy. After verification, the precipitation and total solar radiation show a significant negative mean correlation during the period, with r = −0.64 (p < 0.01). Therefore, the increase of precipitation in the growing season does not promote the growth of vegetation, but instead weakens the photosynthesis by cutting down the energy source from solar radiation.

On the other hand, precipitation shows a positive influence on the variation of NPP in the dry season. This is the result of the extreme lack of water for vegetation growth in these months, when the mean precipitation is 42.58 mm per month, less than one-fifth the level of the growing season. Thus, precipitation replaces radiation as the main limitation for photosynthesis, as increased precipitation could certainly promote vegetation growth by providing the necessary element of water.

4.3.3. Lagged Impact of Precipitation on NPP

The delayed effect of precipitation on NPP should also be considered, and time lags of 0–3 months were considered according to previous studies [82]. From January to December, the correlations between monthly NPP and precipitation in the same month or 1–3 months previous were calculated and compared, as shown in Figure 10. Although NPP in the growing season shows a significant negative correlation with precipitation, optimal positive correlations (meaning the highest correlation coefficients) can be observed with the precipitation of one or two months before. Except for August, during which the optimal time lag is two months, the NPP of the other months in the growing season all show one month lag with precipitation. This indicates that increased precipitation in the previous one or two months can promote the growth of vegetation, because soil can store water from rainfall and provide water to the vegetation in the following months. Meanwhile, this phenomenon cannot be found for the dry season months, during which the optimal correlation between NPP and precipitation appears in the identical month, without lag. This might be the result of the urgent demand for water in these dry months, so vegetation immediately responds to the precipitation, without time lag. This indicates that precipitation has an important impact on the variation of NPP in the study area, even though a negative correlation is observed in the growing season.

4.4. Driving Forces for the Inter-Annual Variation of NPP in the Three Stages

Although precipitation presented three completely parallel stages to NPP, it was not enough to explain the three-stage inter-annual variation of NPP. The variation of NPP in the growing season and dry season was therefore calculated to further study the driving forces for NPP change in the three stages, as shown in Figure 11. It can be seen that the changes of annual NPP in the three periods were dominated by the amounts in the different seasons. For the stage from 1982 to 1992, the decline of annual NPP was apparently dominated by the NPP change in the growing season (slope = −3.34), and the NPP increment from 1992 to 2002 was caused by the NPP variation in the dry season (slope = 4.67). Meanwhile, the decrease of NPP from 2002 to 2014 was the result of the decreased NPP in both the growing season (slope = −1.39) and the dry season (slope = −0.68). The variation trends of the climatic factors in the related seasons and stages were calculated to find the reasons for the NPP change, as shown in

Table 4. Variation trends of NPP and climatic factors in the related seasons and stages.

Time Period	Season	NPP (gC m−2 Year−1)	Precipitation (mm Year−1)	Temperature (°C Year−1)
1982–1992	Growing season	−3.34	−8.25	−0.014
1992–2002	Dry season	4.67	−1.52	0.15
2002–2014	Growing season	−1.39	−9.83	0.054
	Dry season	−0.68	−7.43	0.021

. For the time periods of the growing season from 1982 to 1992 and the two seasons from 2002 to 2014, precipitation showed the same decreasing variation trends as NPP (−8.25, −9.83, −7.43 mm per year, respectively). Precipitation did not show an obvious decreasing trend in the dry season from 1992 to 2002 (only −1.52 mm per year) when the NPP increased. Instead, temperature presented a very significant growth of 0.15 °C per year in this period, which was much greater than the temperature rise in the other stages.

The variation trends of NPP and climatic factors at the pixel level were also calculated and analyzed for the four groups of data. For the area where NPP showed the same variation trends as the mean seasonal NPP, statistics of the climate variation trends were calculated and are shown in

Table 5. The percentage of decreased and increased pixels of precipitation and temperature within the area showing the same variation trends as the seasonal NPP.

Time Period	Season	NPP Trend	Decreased Pixels (%)		Increased Pixels (%)
			P	T	P	T
1982–1992	Growing	Decreasing	78.63	72.66	21.37	27.34
1992–2002	Dry	Increasing	58.92	0	41.08	100
2002–2014	Growing	Decreasing	63.41	5.95	36.59	94.05
	Dry	Decreasing	99.07	26.98	0.93	73.02

Note: P denotes the precipitation; T denotes the temperature.

. Taking the example of the growing season from 1982 to 1992, the percentage of area with increased or decreased precipitation/temperature was counted among the pixels with decreased NPP. For the growing season of the period from 1982 to 1992, 95.54% of the study area showed a decreased NPP in total. Among these regions, 78.63% showed decreased precipitation, with a mean trend of −12 mm per year. The continuous reduction in precipitation was likely responsible for the reduction of NPP, because it would cause persistent droughts. The growth of vegetation was also affected by the continuous decline of temperature in 72.66% of the area, due to the decreased biological activity. The decrease of precipitation would also be the reason for the decline of NPP from 2002 to 2014. Especially in the dry season, 99.07% of the area with a decreased NPP presented a precipitation decrease, with a mean trend of −7.80 mm season⁻¹. The value was 19.21% of the monthly precipitation amount in that period. The percentage of area and mean trend in the growing season also reached 63.41% and −12.22 mm per year. The significant precipitation decreases in each season have caused frequent droughts during this period, which were the reason for the NPP decline, most notably the four-year extreme drought from 2009 to 2012. Meanwhile, for the dry season from 1992 to 2002, although the mean seasonal precipitation slightly decreased, more rainfall occurred in 41.08% of the area. As a result, the rapidly increasing temperature trend of 0.15 °C per year greatly promoted the vegetation growth. The climate warming benefited the vegetation even more in the dry season, when the temperature was relatively low, with a mean value of only 11.03 °C. In general, it can be concluded that the increment of NPP from 1992 to 2002 was mainly caused by the growth of the dry season due to the significant climate warming, and the decline in NPP from 1982 to 1992 and from 2002 to 2014 was due to the frequent droughts caused by the precipitation decrease during these periods.

5. Discussion

NDVI is one of the core parameters in modeling NPP at regional or larger scales, which decides the characteristics of the NPP result. Due to the inter-inhibitive characteristics of current sensors, none of the previous researches could conduct a NPP study with ideal time span and spatial resolution at the same time. Thus, we proposed a new framework to obtain a long-term and suitable resolution NPP series, through fusing multi-source remote sensing NDVI datasets. In this section, we would like to discuss the uncertainties of the fusion results in estimating NPP comparing to the original NDVI data, as well as its applicability to other regions.

5.1. Uncertainties of the Fused NDVI in Estimating NPP

Since the fusion process is the core to estimate NPP, it is necessary to discuss the uncertainties of NPP caused by the fused NDVI. Therefore, the NPP calculated by the fused NDVI (NPP_fused) was compared with the NPP calculated by the original MODIS (NPP_MODIS) and GIMMS3g NDVI (NPP_GIMMS3g), respectively. First of all, the comparison of NPP_fused and NPP_MODIS was conducted in 2011, when the fused NDVI showed a relatively worse result in Figure 4. The result indicated that the NPP_fused shows strong consistency with NPP_MODIS, both in the spatial distribution and the intra-annual variation. The spatial correlation coefficient of the two data reaches 0.94, and the mean absolute difference (MAD) is only 4.16 gC m⁻² month⁻¹. Moreover, their temporal variations are also almost the same, with r as high as 0.99. It means that the fused NDVI has the almost same performances with the true MODIS data in estimating NPP. It is a result of the good data consistency showed in Figure 4 and Figure 5. The two kinds of NDVI data essentially showed little differences, after the multiple operations of CASA model, the NPP results would be even more consistent.

Meanwhile, in order to demonstrate the impacts of the NDVI data processes in NPP estimation, the NPP based on the original GIMMS3g NDVI (NPP_GIMMS3g) and the MODIS-like NDVI (NPP_like) were also calculated and validated with the in situ data in 1983. As shown in Figure 12, and compared to the performances of NPP_fused in Figure 6, the NPP_fused is obviously more accurate than the other two results. The NPP_GIMMS3g is the worst one, with overestimation in almost all the plots (r = 0.33). The performance of the NPP_like is better (r = 0.44), due to the NDVI being normalized to be consistent with the MODIS data, which is more accurate with lower values. The promoted spatial resolution would be the reason for the better accuracy of NPP_fused, because it could reduce the overestimation and precision loss caused by the spatial heterogeneity when modeling NPP.

In general, according to the evaluations above, the fused NDVI would not bring extra uncertainties in estimating NPP. The fused NDVI showed the extremely similar performances with true MODIS data in estimating NPP, and could obtain a much more accurate NPP than the original GIMMS3g NDVI data.

5.2. Applicability of the Fusion Method to Other Regions

This study only applied the new framework to estimate NPP in Yunnan province, to investigate the regional NPP patterns and its relation to the complex climate. Thus, it is necessary to discuss the applicability of the new framework to estimate NPP in other regions. The core of the framework is the generation of a new long-term NDVI time series with higher spatial resolution, by fusing multi-source remote sensing data. Thus, we can prove the feasibility of our framework, if the fusion processes are applicable in other regions. In this case, we additionally selected the typical and important terrestrial ecosystem in Qinghai-Tibet Plateau and Amazon area, as the testing region, and the same NDVI processes in Section 3.2 was conducted. The tile numbers of the MODIS data are ‘h11-v09’ and ‘h24-v05’, respectively. After the fusion processes, the same simulated validation in the Section 4.1.1 was applied. The quantitative evaluation is shown in Figure 13, and we could observe the same conclusions as in the study area of Yunnan. The fused NDVI presents good consistency with the true MODIS data, and the fine results maintain a stable accuracy as the years between the prediction date and reference date increase. The value of the statistics shows small differences, which would be the result of the different levels of NDVI value in different regions. Besides, also similar to the conclusions in Yunnan, the qualitative comparisons indicate a clearly improved spatial resolution for the fused NDVI, and showing a highly consistent spatial distribution with the true MODIS data. From the experiments above, it can be concluded that the fusion processes could obtain fine results in different regions with different surface conditions. Therefore, the same framework described in this paper should be applicable in other regions, to obtain the long-term NPP series with higher spatial resolution.

In summary, this study proposed an efficient data-processing framework to map a long-term NPP with enhanced spatial resolution, at regional scale. The method would not add extra uncertainties, but could improve the NPP accuracy by weakening the error caused by spatial heterogeneity. It agrees to the previous studies, which have proved that coarse resolution can lead to an obvious overestimation and accuracy loss when modeling NPP [38,39,40]. Furthermore, the three stages of the NPP annual variation in this study also agree to the result of Liang et al. [83], some differences in break points and trends value would be the result of the different study regions. However, further efforts should be made to investigate the more detailed spatial and temporal patterns of NPP in such complex terrain, such as the specific variation trends and the response of NPP to climate in different elevation zones. Furthermore, the relationships between NPP and more environmental factors (i.e., land surface parameters, other climatic factors, atmospheric condition, land-use/land-cover changes, disaster disturbance, human activities, and so on) also need to be investigated, comprehensively considering the linear and non-linear relationships or trends, in order to obtain a better understanding of the influencing mechanism of the regional carbon cycle.

6. Conclusions

In this paper, we conducted a study on the regional carbon cycle in Yunnan province of southwest China, from 1982 to 2014, at a 1-km scale. As a core parameter for NPP estimation, the NDVI time series were innovatively obtained by fusing the GIMMS3g data and MODIS data to combine their respective advantages. The comprehensive evaluations demonstrated that the fused NDVI is as accurate as the MODIS NDVI for NPP estimation, and could greatly improve the NPP accuracy compared to the original AVHRR NDVI. The final estimated NPP presented good consistency with the field measurements, both in all (r = 0.79) and in each vegetation types.

We analyzed the NPP series and got some interesting conclusions. NPP generally decreased from southwest to northeast in spatial terms, and more than 67% of the study area showed an increasing NPP in the last 33 years. The Annual NPP variation could be separated to three distinct stages, and its overall increasing trend was mainly contributed to by the growth during 1992–2002. There was a significantly negative correlation between precipitation and NPP in the growing season, but one or two months lagged positive correlations were observed. However, in the dry season, the correlations were positive without time lag. This huge seasonal heterogeneity led to a non-negative correlation at the annual scale, because NPP in the growing season it was far greater than in dry season. Furthermore, the frequent droughts caused by the precipitation decrease led to the NPP decline during 1982–1992 and 2002–2014; and NPP rapidly increased from 1992 to 2002 due to the significant climate warming, when the precipitation varied only slightly. Finally, we discussed the uncertainties and the transferability of the fusion framework in estimating NPP; the result indicated that the framework could improve NPP accuracy by enhancing the spatial resolution and could be efficiently applied in other regions.