The Temporal Analysis of Regional Cultivated Land Productivity with GPP Based on 2000–2018 MODIS Data

Abstract

The spatiotemporal change characteristics of Cultivated Land Productivity (CLP) are imperative for ensuring regional food security, especially given recent global warming, social development and population growth. Based on the hypothesis that the Gross Primary Productivity (GPP) is a proxy of land productivity, the Moderate Resolution Imaging Spectroradiometer (MODIS) data with 500-m spatial resolution and 8-day temporal resolution was employed by the Vegetation Photosynthesis Model (VPM) to calculate GPP in Jilin Province, China. We explored the level of CLP using the GPP mean from 2000 to 2018, and analyzed the changing trend and amplitude of CLP in the whole study period using both Theil–Sen median trend analysis and the Mann–Kendall (MK) test, and forecasted the sustainability of CLP with the Hurst exponent. The trend result and the Hurst exponent were integrated to acquire the future direction of change. The results revealed that: (1) The CLP level was generally high in the southeast and low in the northwest in cultivated land in Jilin, China. The area with the lowest productivity, located in the northwest of Jilin, accounted for 15.56%. (2) The majority (84.77%) of the area showed an increasing trend in 2000–2018, which was larger than the area that was decreasing, which accounted for 3.97%. (3) The overall change amplitude was dominated by a slightly increasing trend, which accounted for 51.48%. (4) The area with sustainability accounted for 33.45% and was mainly distributed in the northwest of Jilin. The area with anti-sustainability accounted for 26.78% and was mainly distributed in the northwest and central Jilin. (5) The Hurst exponent result showed that uncertain variation of CLP is likely to occur in the future over the entire region, and the central region is prone to display degeneration. Therefore, the results of this study indicated that quality improvement policy could be implemented for the middle-to-low yield fields in northwest Jilin, and dynamic monitoring and protection measures could be implemented for the areas with uncertain future changes and decreasing sustainability.

1. Introduction

The CLP plays a critical and irreplaceable role in maintaining regional sustainable development and ensuring national food security [1,2,3]. With increasing food demand, it is urgent to explore CLP dynamics spatially and temporally at national and regional scales [4]. However, many studies have been conducted to research the spatial distribution of CLP [5,6], such as land productivity potential assessment [7] and soil quality evaluation [8]. In the context of increased soil erosion and climate change [9], compared with spatial assessment studies, understanding the temporal trend and future sustainability of CLP has become the key to regional sustainable development [10,11]. It is also essential for improving the level of CLP and developing cultivated land protection policies to control land and environmental degradation [12]. The research and development of temporal changes of CLP are limited, mainly for the follow reasons. The first is that it is difficult to find an indicator of CLP which characterizes the combined effect of various factors over many years [13]. The other reason is that the trend and amplitude of temporal change are difficult to quantify. Furthermore, how to determine future sustainability is also a problem [14]. Therefore, some related research had been carried out to address these problems.

Cultivated land is a complex system that combines natural and human factors [15]. At present, studies on the indicators for monitoring and estimating CLP include the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), GPP which reflects crop growth based on remote sensing (RS) [16,17,18]. These indicators have a high correlation with crop yield. As far as the comprehensiveness of these indicators is concerned, the effect of temperature and precipitation on crop growth are considered in the GPP calculation, which has a higher accuracy of simulated production [19]. The GPP is the capacity of the vegetation to capture carbon and energy by photosynthesis over a unit of time [20]. It is a comprehensive indicator formed by the interaction of various factors, such as climate, topography, soil, plant and microbial characteristics, disturbances and human influences. Hence, GPP is considered one of the major indicators in agricultural ecosystems for monitoring and assessing the variations in productivity [21,22]. The models that estimate the GPP of terrestrial ecosystems have been grouped into three categories [23]: process-based models, regression statistical models and satellite-based drive models. The advantage of the satellite-based drive model compared to other models is that it uses simple parameters and large-scale inversion. RS technology, with long temporal series and high spatial coverage, has developed into an indispensable and useful method for monitoring land surface dynamics and changes [24]. VPM, a satellite-based drive model of GPP estimation, is developed from RS data and flux observations, with the advantage of simple input parameters and equal accuracy with process-based models. The most important merit is that the C₃ and C₄ crops of different photosynthesis processes are considered separately. Since 2004, model validation studies have been carried out in 10 types of ecosystems at 21 field sites [25]. The VPM was used as an important basis for GPP time trend extraction based on RS data [26]. The studies of monitoring and analyzing trend changes with long-term RS data have mainly focused on the dynamic monitoring of inter-annual vegetation [27]. Time-series statistical methods, such as Theil–Sen median trend analysis and MK test, T-test, linear regression analysis, Cusum test, Pearson correlation coefficient, Spearman’s rho test, Wald–Wolfowitz test and the Hurst exponent [28], were used to extract dynamic information. Some hypothesis conditions are required in some statistical methods of data, such as the normal distribution of data or a large number of sample data, and the statistical results are subject to outlier interference [29]. Theil–Sen median trend analysis and the MK test, which are simple in calculation, do not need to be submitted to specific statistical distribution, and are not affected by outliers, are widely used in time series trend analysis of runoff, temperature, precipitation and water quality. The combination of the two methods can not only reflect the changing trend but also quantify the amplitude of change. The Hurst exponent [30] has been used to judge whether the time series data belong to a random walk or a biased random walk, which has been used to analyze the sustainability characteristics of vegetation growth trends [31].

Jilin Province is the largest commodity grain production base in China and plays a vital role in ensuring Chinese food security. Since the implementation of Chinese commercial grain base construction policy from 1983 to 2004, 33 commercial grain base counties have provided 276 billion km of grain for China, accounting for 10% to 15% of the total Chinese commercial grain [32]. With the increase in population demand for grain, the level of CLP in Jilin, which is the main granary in northeastern China, has been a focus. However, there are increasing rates of change in climate and deteriorating environmental conditions in Jilin. Over the past 50 years, the average temperature has risen by 0.35 ℃/decade and the average precipitation has dropped by 16.34 mm/decade in this region [33]. In the last decade years, the organic matter content (OMC) of black soil decreased from 6% to 2% in most of Jilin, which has caused some black soil farmland to lose production capacity [34]. Fertilizer is being applied 10 times more than 30 years ago and the utilization of pesticides has increased from 2.2 tons to 25,000 tons, which is the main reason for soil compaction and soil pollution [35]. It is necessary to determine the changing trend and variation amplitude of CLP in Jilin under the background of climate change and soil degradation. Further, the prediction of changes in the short term is the key to optimize the allocation of cultivated land resources, to maintain the stability of productivity and enhance the low- and medium-yield fields in Jilin.

This paper is based on the hypothesis that the GPP is a proxy of land productivity, and in this paper, we analyze the CLP level, the trend and amplitude of change in Jilin from 2000 to 2018 and the future direction of CLP. The GPP time series is estimated by the VPM model derived from MODIS data. The level of CLP was estimated using the GPP mean from 2000 to 2018. The trend and amplitude of GPP were quantified by both the Theil–Sen median trend analysis and the MK test. The sustainability characteristic of future change was assessed by the Hurst exponent. The MK trend and Hurst exponent were integrated to express the future direction of change. The objectives were to determine: (1) the temporal and spatial characteristics of the CLP level; (2) the trend and amplitude of CLP during 2000–2018; and (3) the sustainability and variation direction of GPP in the short-term future. Based on this study, the level of regional CLP improvement and cultivated construction achievements are promoted by understanding the temporal variation and the spatial distribution of CLP in the Jilin.

2. Materials and Methods

2.1. Study Area

The Jilin Province is located in the middle of northeast China, from 121°38′ to 131°19′ E and 40°50′ to 46°19′ N, with a total area of 18.74 × 10⁴ km², as shown in Figure 1. Jilin Province has a temperate continental monsoon climate with four distinct seasons, with rain and heat in the same season. The study area is divided into the eastern mountainous area and the western plain by the central Black Mountain. The cultivated land is about 5.54 million ha, accounting for 28.89% of the Jilin Province area, 70% of which is distributed in the central and western plains. The cultivated land in the Jilin is mainly dry land and paddy fields. According to the difference in the geographical environment, the study area was divided into the northwest ecologically fragile area, the central area with high quality and contiguous cultivated land and the southeast area located in the mountainous area with scattered distribution. The main crops of Jilin Province are C₄-based maize (Zea mays L.) and C₃-based paddy rice (Oryza sativa L.).

2.2. Data Source and Pre-Processing

The 8-day composite MOD09A1 dataset of MODIS was downloaded from the USGS website (www.usgs.gov/). The MOD09A1 is a level 3 product from the Terra platform with 500 m spatial resolution, where geometrical corrections and cloud contamination removal were done. The quality assessment band-flagged pixels of high quality were used for the analysis. A total of 380 images from 1 May 2000 to 1 October 2018 were applied for GPP estimation. The HDF format was transformed into GeoTIFF format and the sinusoidal projection was converted to the WGS84 projection with the MODIS Reprojection Tool (MRT), which can be acquired from the NASA Land Processes Distributed Active Archive Center. The cultivated land grids with 500 m spatial resolution were collected from the agricultural land classification database. The images were masked by the cultivated land grids using the Google Earth Engine (GEE) platform. The meteorological datasets consisting of 8-day mean temperature data were collected from the National Meteorological Center of China. The data were from 51 meteorological stations between 2000 and 2018. The observations were interpolated over 8 days with the same resolution as the MOD09A1 for obtaining spatiotemporal distributions of temperature to estimate the GPP. The flowchart of this paper is shown in Figure 2.

2.3. Methods

To provide temporal characteristics of the CLP in the Jilin Province, the following investigations were conducted:

• The VPM model was used to estimate GPP for C₃-based paddy rice and C₄-based maize. Then, two types of comparable GPP were obtained using an anomaly analysis and the GPP time series was generated.
• To investigate the temporal and spatial characteristics of CLP over a long period, the level of CLP was analyzed by GPP multi-years mean. The trend and amplitude were quantified by both the Theil–Sen median trend analysis and the MK test.
• The sustainability characteristic of future change was assessed by the Hurst exponent, which was integrated with the trend result to express the future direction of change.

2.3.1. GPP Estimation with the VPM Model and the GPP Mean Calculation

The VPM model is a GPP estimation model driven by the RS indexes, such as the EVI and land surface water index (LSWI), which was developed by Xiao [36] and is based on conceptual partitioning of vegetation canopies into nonphotosynthetic vegetation (NPV) and photosynthetically active vegetation (PAV) or chlorophyll (CHL). The EVI and LSWI were calculated by the MOD09A1 bands that include blue (459–479 nm), red (620–670 nm), near-infrared (841–875 nm) and shortwave infrared (SWIR, 1628–1652 nm). The calculation formulas of the EVI and LSWI and the calculation procedures of the model are referred to in [23]. Two indexes were combined with other datasets as input parameters of the VPM model to calculate the GPP time series data (Figure 3.). The formula is as follows.

$$GPP={\epsilon }_0\times W_{scalar}\times P_{scalar}\times T_{scalar}\times FAPA{R}_{PAV}\times PAR$$

(1)

where GPP is the cultivated land GPP (Kg·C/m²), ${\mathsf{\epsilon }}_0$ is the maximum light use efficiency, and T_scalar, W_scalar, P_scalar are the scalars for the effects of temperature, water, and leaf phenology on light use efficiency of vegetation, respectively. PAR (MJ/m²) is the photosynthetically active radiation in a period, FAPAR_PAV represents the fraction of PAR absorbed by the vegetation canopy.

To obtain the spatial distribution of the level of CLP during 2000–2018, the GPP mean of the period was calculated by the following formula:

$$\overline{GPP}=\frac{1}{n}{\displaystyle \sum }_{i=1}^{n}GPP$$

(2)

where $\overline{GPP}$ is the mean during 2000–2018, and n is the number of years.

2.3.2. Theil–Sen Median

The Theil–Sen median was used to extract the trends of the GPP time series. Theil–Sen median trend analysis is a robust trend statistical method that uses the median slopes between all n(n − 1)/2 pair-wise combinations of the time series data, which contribute to reducing the influence of missing or abnormal data and is particularly effective for the estimation of the trend in a time series [37]. The slope of the Theil–Sen median represents an increase or decrease with the GPP series data between 2000 and 2018 on a pixel scale. The calculation of the slope was as follows:

$$Slop{e}_{GPP}=Median\left(\frac{GP{P}_j-GP{P}_i}{j-i}\right),2000\le i\le j\le 2018$$

(3)

where i, j are the number of years in the study period, and GPP_i and GPP_j represent the GPP values at times i and j, respectively. An increasing trend is inferred from $Slop{e}_{GPP}$ > 0, and a decreasing trend is inferred from $Slop{e}_{GPP}<0$. The significance of the trend was quantified using the MK test.

2.3.3. The MK Test

The MK test, which is a non-parametric statistical method that does not require compliance with a certain distribution and is not disrupted by other abnormal values, was used to quantify the significance of the trend. In the MK test, the null hypothesis (H₀) explains that the sequence (GPP₂₀₀₀, GPP₂₀₀₁, GPP₂₀₁₈) consisted of a sample of n independent and uniformly distributed random variables. The alternative hypothesis (H₁) of the two-sided test was that the distribution of x_j and x_k is not identical for all values of j and k. S was calculated as follows:

$$S={\displaystyle \sum }_{k=1}^{n-1}{\displaystyle \sum }_{j=k+1}^{n}Sgn\left(GP{P}_j-GP{P}_k\right)$$

(4)

where S is the test statistic with a normal distribution. GPP_k are ranked from k = 1, 2, …, n – 1 and GPP_j from j = k + 1, …, n. Each data point of GPP_k is used as a reference point and is compared with all other data points of GPP_j:

$$\mathrm{Sgn}\left(GP{P}_j-GP{P}_i\right)=\left[\begin{array}{c}1\left(GP{P}_j-GP{P}_i\right)>0\\ 0\left(GP{P}_j-GP{P}_i\right)=0\\ -1\left(GP{P}_j-GP{P}_i\right)<0\end{array}\right]$$

(5)

$$\mathrm{Z}=\{\begin{array}{c}\frac{S-1}{\sqrt{VAR\left(S\right)}},S>0\\ 0,S=0\\ \frac{S+1}{\sqrt{VAR\left(S\right)}},S<0\end{array}$$

(6)

$$\mathrm{VAR}(\mathrm{S})=(n(n-1)(2n+5)-{{\displaystyle \sum }}_{i=1}^m{t}_i(t_i-1)(2t_i+5))/18$$

where the range of Z is (−∞, +∞), and $\left|\mathrm{z}\right|>{\mathsf{\mu }}_{1-\mathsf{\alpha }/2}$ is the significance level, while the sign of this time series shows significant changes on the level of $\mathsf{\alpha }$. In this study, the $\mathsf{\alpha }$ was set as 0.05 to analyze the significance of the GPP trend from 2000 to 2018 at a confidence level of 0.05.

2.3.4. The Hurst Exponent

The Hurst exponent was proposed by British hydrologist Hurst (1951), and later revised by Mandelbrot and Wallis (1969), for determining the sustainability of the time series in the future. This method was used to assess the sustainability of vegetation growth by detecting whether the NDVI has long-term dependence [37,38]. In this paper, the sustainability of CLP in the short term was judged by the quantitative description of GPP dependence using the Hurst exponent. The calculation steps are as follows:

(1)

To define the time series GPP(t), t = 1, 2...n.

(2)

To calculate the mean of the GPP time series,

$${\overline{GPP}}_{\left(\tau \right)}=\frac{1}{\tau }{{\displaystyle \sum }}_{t=1}^{\tau }GP{P}_{\left(t\right)}\tau =1,2,\dots ,n.$$

(7)

(3)

To calculate the accumulated deviation,

$${\mathrm{X}}_{\left(t,\tau \right)}={{\displaystyle \sum }}_{t=1}^t\left(GP{P}_{\left(t\right)}-{\overline{GPP}}_{\left(\tau \right)}\right)1\le t\le \tau .$$

(8)

(4)

To acquire the level difference,

$${\mathrm{R}}_{\left(\tau \right)}=\underset{1\le t\le \tau }{\max }GP{P}_{\left(t,\tau \right)}-\underset{1\le t\le \tau }{\min }GP{P}_{\left(t,\tau \right)}\tau =1,2,\dots ,n.$$

(9)

(5)

To acquire the standard deviation sequence,

$${\mathrm{S}}_{\left(\tau \right)}={\left[\frac{1}{\tau }{{\displaystyle \sum }}_{t=1}^{\tau }{\left(GP{P}_{\left(t\right)}-GP{P}_{\left(\tau \right)}\right)}^2\right]}^{1/2}\tau =1,2,\dots ,n.$$

(10)

(6)

To acquire the H exponent,

$$\frac{R\left(\tau \right)}{S\left(\tau \right)}={\left(c\tau \right)}^{\mathrm{H}}.$$

(11)

According to previous research [39], the value of H ranges from 0 to 1, and there are three types: When H = 0.5, it indicates that the GPP time series is a random sequence without sustainability, and it is inferred that the trend of the CLP in the future is unstable. When H > 0.5, it means that the GPP time series is a sequence with sustainability, the future trend is consistent with past changes, and the trend of the CLP in the future is sustainability. When H < 0.5, it means that the GPP time series is a sequence with anti-sustainability, the future trend is opposite to past changes, and the trend of the CLP in the future is anti-sustainability.

3. Results

3.1. The Characteristics of the Temporal and Spatial Distribution of the CLP Level

The level of CLP was accessed by removing the inter-annual variation difference of GPP from 2000 to 2018 (shown in Figure 4a,b). According to the classification of the natural breakpoint (

Table 1. CLP classification standard and area statistics.

GPP (Kg·C/m2)	CLP Grade	Area Percent (%)
0~0.45	Level 1	15.56
0.45~0.6	Level 2	6.37
0.6~0.7	Level 3	38.43
0.7~0.88	Level 4	30.84
0.88~1.29	Level 5	8.80

), the CLP level was classified into levels 1–5, from low to high. The whole was 0.633 Kg·C/m², which means it was level 3, called the provincial level. The temporal trend of the annual GPP during 2000–2018 is shown in Figure 5a. The whole area fluctuated in the range of 0.57–0.74 Kg·C/m² and exhibited an upward trend with a slope of 0.0063 Kg·C/m²yr in the monitoring periods. The maximum was in 2017 with 0.57 Kg·C/m² and the minimum was in 2000 with 0.74 Kg·C/m². Temporal periodicity was seen with a five-year cycle of increasing and then decreasing. More specifically, the temporal change at different levels is shown in Figure 5b, where it is evident that all slopes were higher than 0, which means that they all showed an increasing trend from 2000 to 2018. However, the slope rate of each level was different, and the fastest change rate was 0.0086 Kg·C/m²yr in level 5. The slowest change rate was with an annual growth rate of 0.00427 Kg·C/m² in level 1. Other levels fluctuated within this range during this period. The important conclusion was drawn that a higher level of productivity had a slower rate of increase and the level with the lower productivity had a faster rate of increase.

The spatial distribution of the CLP level is shown in Figure 4 and Figure 5c. The spatial distribution characteristics of the northwest at low altitude and southeast at high altitude were displayed in the cultivated land in Jilin. The level 3 area was the most abundant and accounted for 38.4% (Table 1). It was distributed in the central region, such as Siping City and Changchun City. The level 2 area accounted for 6.37%, which was the lowest, and was concentrated in Baicheng City and Songyuan City in the northwest region. The level 1 area occupied 15.56% of the total area, and was mainly concentrated in the northwestern part of Jilin, such as Baicheng City and Songyuan City. The level 4 area occupied 30.84% of the total area, and was mainly concentrated in the southeastern part of Jilin, such as Jilin City, Liaoyuan City, Tonghua City and Yanbian City. The level 5 area accounted for 8.8% of the total area. This area was mainly found alongside the level 4 area, including southeast Jilin City, north Tonghua City, Baishan City and Yanbian City.

3.2. The Trend and Amplitude of Change of the GPP from 2000 to 2018

Theil–Sen median trend analysis and the MK test were used to detect the trend and amplitude for each pixel. The results were divided into three categories according to the actual situation of the S_GPP. There was no area where the S_GPP was 0. Where the S_GPP was between −0.0005 and 0.0005 were classified as stable areas. Where the S_GPP was greater than or equal to 0.0005 were classified as increasing areas. Areas with an S_GPP less than −0.0005 were classified as decreasing areas. Moreover, the significance test results of the MK test were used to assess the amplitude of change. At the confidence level of 0.05, the area was classified as significant variations (Z > 1.96 or Z < −1.96) or slight variations (−1.96 ≤ Z ≤ 1.96). By integrating the classification results of the Theil–Sen median trend analysis and the MK test, it was comparable to the data of trend variations of the CLP and was classified into five categories.

Areas percentages of CLP trends and their distributions are shown in

Table 2. Statistics of different trend areas from 2000 to 2018.

SGPP	Z	Trend	Percent (%)
$\ge 0.0005$	$\ge 1.96$	Significant increase	33.29
$\ge 0.0005$	−1.96–1.96	Slight increase	51.48
−0.0005–0.0005	−1.96–1.96	Stable	11.26
$<-0.0005$	−1.96–1.96	Slight decrease	2.51
$<-0.0005$	<−1.96	Significant decrease	1.46

and Figure 6, respectively. (1) The areas where the CLP showed a significantly increasing trend occupied 33.29% of the total area. They were mainly distributed in Baicheng City and Songyuan City in the northwest region, but were also scattered in the west Siping City and northwest edges of Changchun City in the central region. (2) The areas where the CLP showed a slight increasing trend occupied 51.48% of the total area. These areas were distributed throughout the region and mainly concentrated in Changchun City, central Siping City, scattered in northwest Baicheng City, central Songyuan City in the northwest region and central Liaoyuan City and north Tonghua City in the southeast region. (3) The areas where the CLP showed a stable trend accounted for 11.26% of the total area. These areas were mainly found alongside the areas with a slightly increasing trend, including northwest Yanbian City and southeast Jilin City. (4) The areas where the CLP showed a slightly decreasing and a significantly decreasing trend accounted for 2.51% and 1.46%, respectively. These areas were scattered in Yanbian City, Liaoyuan City, and Tonghua City in the southeast region (Figure 6c). The amplitudes of all types are shown in Figure 7. The rates of significant changes, slight changes and no change were 0.00739 Kg·C/m²·yr, 0.00398 Kg·C/m²·yr, and 0.00113 Kg·C/m²·yr, respectively.

3.3. The Sustainability and Direction of Change of the CLP in the Future

As shown in Figure 8, the average Hurst exponent H value of GPP in Jilin Province was 0.52, which was greater than 0.5 and indicates that the CLP change in the future will be consistent with the past. Namely, it will be sustainable. As far as the spatial distribution of the partition was concerned, the area with sustainability accounted for 33.45% and was mainly distributed in Baicheng City and Songyuan City in the northwest, and Liaoyuan City in the southeast, of Jilin province. An H value less than 0.5 suggested that the future direction of change of the CLP will be opposite to the past, namely, anti-sustainability, which accounted for 26.78% and was mainly distributed in the northwest of Baicheng City, Siping City and Changchun City in the central Jilin. The H value equal to 0.5 suggested that the future direction of change of the CLP will be random, namely, stable, which accounted for 60.23% and was mainly distributed in the northwest of Changchun City, Liaoyuan City and Jilin City in Southeast Jilin. Overall, the stable and anti-sustainability areas showed signs of large fluctuations in CLP that need to be monitored.

To reveal the direction of change of CLP in the future, the Theil–Sen S value, which characterized the past trend, was superimposed by the H value and characterized the sustainability of future changes to generate the comprehensive information of the trend and sustainability. The results were summarized into six categories from ‘sustainability increased’ to ‘uncertain’ as shown in

Table 3. Statistical result of the trend and the Hurst exponent.

SGPP	H	Variation Types	Percentage (%)
>0.0005	>0.5	sustainability increased	13.37
<−0.0005	>0.5	sustainability decreased	13.98
>0.0005	<0.5	Anti-sustainability increased	14.54
<−0.0005	<0.5	Anti-sustainability decreased	12.24
−0.0005–0.0005	>0.5	Stable	6.1
----	=0.5	uncertain	39.77

and Figure 9. The area shown to be uncertain in the future was the largest in the Jilin Province and accounted for 39.77%, which was mainly distributed in the Changchun City and Siping City in the central part. The ‘sustainability increased’ and ‘anti-sustainability decreased’ areas, which indicate that CLP will be improved in the short term, accounted for 25.61% and were mainly located in Baicheng City and Songyuan City. The ‘sustainability decreasing’ and ‘anti-sustainability increasing’ regions, where the CLP deteriorate in the short term, accounted for 28.52% and were mainly located in the southeastern City and Songyuan City in the central part. The area that was shown to be stable in the future accounted for 6.1%, and was scattered in the north of Baicheng City.

4. Discussion

In terms of different regions, the level of CLP, the amplitude of the trend and the future sustainability were combined with the geographical characteristic to explore the driving forces that are contributing to guiding the CLP improvement projects. In the northwestern region of Jilin Province, the CLP level is low at level 2 (Figure 5b) and a significant increase was shown during 2000–2018, which will keep sustainability increasing in the future. The average annual precipitation in this area is 400–700 mm, which is less than the 700–900 mm in the southeast of Jilin. Most of the cultivated land was low-lying, which will degenerate to saline-alkali land and wind-blown land, where the ecological environment is fragile. It is also more vulnerable to natural disasters such as flooding, salinization and degradation, and the artificial irrigation facilities are poor in this area. These factors limited the regional CLP level. As shown in Figure 10, precipitation has greater regional differences compared to temperature in the northwest, than in the central part. Therefore, precipitation is more likely to affect regional productivity levels than temperature. In recent years, measures such as the construction of a cultivated land shelterbelt, high-standard cultivated land construction and water-saving irrigation projects have been implemented in the central and western regions to effectively control disasters and improve environmental conditions, which might be the reason for the significant increase in the region. The improvement of climatic conditions such as the increase of precipitation and the rise of temperature over the years is conducive to agricultural production, while the development of cultivated land regulation projects has gradually improved the water and fertilizer conservation capacity of cultivated land. At the same time, due to the continuous improvement in the technical production level, such as the cultivation of drought-resistant and high-yield seeds and complete cultivated land management technology, the CLP level in the northwest region presents a sustainable increasing trend in the future. In the central region of Jilin Province, the average annual productivity is level 3 and it showed a slight increase during 2000–2018, whereas unsustainability decreased and the future was found to be uncertain. The precipitation in the central region is higher than that in the northwest (shown in Figure 10). The artificial irrigation facilities are sufficient to meet the regional needs. Therefore, the productivity level has been higher in the central region for about twenty years. With the increase in precipitation and temperature in recent years [40], productivity is gradually increasing slightly. The average annual productivity level in the southeastern part of Jilin is at level 4. There were a slight increase and stable trend for 2000–2018, and a decrease in sustainability in the southeast region. The cultivated land in the southeast is mainly scattered sloping land. The soil humus is rich and the OMC is high. The precipitation is the highest in the whole region, and so the GPP is the largest. As the area with the highest CLP level in the whole region, with the increase of precipitation and temperature in the future, it may exceed the suitable temperature for crop growth. At the same time, a broken field in the south makes it less resistant to disasters, which may lead to a decrease of CLP in the future.

The 500-m MODIS data can give very accurate surface information. Moreover, the cultivated land in Jilin Province is large and concentrated, so there are less mixed pixels in the 500-m MODIS data. Therefore, it is reasonable to use MODIS to study cultivated land in Jilin. The GPP was selected as the indicator of CLP because, compared with NDVI, the estimated model is combined with water, temperature, and photosynthetic radiation effective coefficient indicators, which are the conditions inherent in the productivity of cultivated land. The method combining the Theil–Sen median slope analysis and the MK test is superior to linear regression analysis, as it can not only obtain trend information but also quantify the amplitude of the trend, and furthermore, it has an early warning effect on preventing the degradation of productivity and appropriately reducing cultivated land-use intensity. In addition, the time series data do not need to obey specific distributions, strong error avoidance, and significance level tests. The sustainability, shown by the H value, only indicates whether the future changes are consistent with the past or not. The direction of change is still unclear. Therefore, the future CLP trend was quantified by superimposing the trend analysis results and the Hurst exponent. The result of coupling these tests is that it reflects the impact of past trends on the future more scientifically.

Based on the GPP time series of Jilin Province from 2000 to 2018, the trend in spatial and temporal variation, the amplitude of change and the sustainability of future changes of CLP were studied. However, different crop types may exhibit different inter-annual variations. The main crop types in Jilin are rice, corn, soybean, sorghum, etc. Therefore, more detailed work is needed to explore the effects of different crop types on GPP temporal and spatial changes. Based on studying the spatial and temporal changes of GPP for many years, combined with multi-source data such as soil type, high-standard farmland construction data, and other social and climate data, we could explore the driving force analysis of the trend, and provide effective regulation for providing regional productivity degradation policy. While the changes in GPP are affected by climate change, human interference, and national or local land and agricultural policies, the drivers associated with the CLP should be explored further. These explorations will be helpful to fully understand the regional CLP, and local governments should adopt policies or projects that promote regional CLP.

5. Conclusions

In this paper, the GPP time series data from 2000 to 2018 were estimated by the VPM based on 500-m MODIS MOD09A1 data. The multi-year mean, the Theil–Sen Median trend analysis, the MK test and the Hurst exponent methods were used to explore the CLP level, the trend and the amplitude of change, and the future sustainability and future direction of change of the regional CLP. The following conclusions were drawn:

(1) Overall, the CLP level in Jilin Province was characterized during 2000 and 2018 by a medium level GPP. The CLP of Jilin was generally high in the southeast and low in the northwest. Therefore, when the cultivated land protection and land improvement measures were implemented, the northwest region with low productivity should be given more consideration.

(2) Most of the regions showed an increasing trend throughout 2000–2018, which was larger than that of the decreasing trend. Most changes were slight. Higher productivity had a lower rate of increase and lower productivity had a higher rate of increase.

(3) The area shown to be uncertain in the future is located in Changchun City and Siping City in the central part. These areas should be monitored dynamically. The improving areas were located in Baicheng City and Songyuan City, and should be given protection measures. The deteriorating area was mainly located in the southeastern cities and the central part, and should undergo comprehensive quality improvement measures.

(4) We propose that the temporal analysis methods are suited for exploring the trend and forecasting the direction of change in the future. Because they are universal, the methods could also be applicable for a finer data source and for other areas in the world to explore future trends.