Detecting the Phenological Threshold to Assess the Grassland Restoration in the Nanling Mountain Area of China

Abstract

The dynamics of vegetation changes and phenology serve as key indicators of interannual changes in vegetation productivity. Monitoring the changes in the Nanling grassland ecosystem using the remote sensing vegetation index is crucial for the rational development, utilization, and protection of these grassland resources. Grasslands in the hilly areas of southern China’s middle and low mountains have a high restoration efficiency due to the favorable combination of water and temperature conditions. However, the dynamic adaptation process of grassland restoration under the combined effects of climate change and human activities remains unclear. The aim of this study was to conduct continuous phenological monitoring of the Nanling grassland ecosystem, and evaluate its seasonal characteristics, trends, and the thresholds for grassland changes. The Normalized Difference Phenology Index (NDPI) values of Nanling Mountains’ grasslands from 2000 to 2021 was calculated using MOD09A1 images from the Google Earth Engine (GEE) platform. The Savitzky–Golay filter and Mann–Kendall test were applied for time series smoothing and trend analysis, and growing seasons were extracted annually using Seasonal Trend Decomposition and LOESS. A segmented regression method was then employed to detect the thresholds for grassland ecosystem restoration based on phenology and grassland cover percentage. The results showed that (1) the NDPI values increased significantly (p < 0.01) across all grassland patches, particularly in the southeast, with a notable rise from 2010 to 2014, and following an eastern to western to central trend mutation sequence. (2) the annual lower and upper NDPI thresholds of the grasslands were 0.005~0.167 and 0.572~0.727, which mainly occurred in January–March and June–September, respectively. (3) Most of the time series in the same periods showed increasing trends, with the growing season length varying from 188 to 247 days. (4) The overall potential productivity of the Nanling grassland improved. (5) The restoration of the mountain grasslands was significantly associated with the grassland coverage and mean NDPI values, with a key threshold identified at a mean NDPI value of 0.5 for 2.1% grassland coverage. This study indicates that to ensure the sustainable development and conservation of grassland ecosystems, targeted management strategies should be implemented, particularly in regions where human factors significantly influence grassland productivity fluctuations.

1. Introduction

The grassland ecosystem in China acts as a vital green ecological barrier, providing essential ecological services such as climate regulation, water conservation, soil improvement, wind prevention, sand consolidation, and biodiversity preservation [1]. Covering approximately 2.64 million km² and accounting for over 27.5% of the country’s land area [2], China’s grasslands can be classified into four distinct ecological functional regions based on their geographical differences, which are driven by differing climates and other natural conditions [3]. Nanling, one of the four ecological regions, has experienced long-term grassland degradation as a result of various natural and anthropogenic factors [4]. This degradation not only impacts the development of animal husbandry and the livelihoods of farmers and herders, but it also poses a serious threat to national ecological security by potentially triggering various natural disasters [5].

Studying vegetation phenology plays a pivotal role in understanding interannual changes within terrestrial ecosystems. It serves as a critical indicator of climate–vegetation interactions, variations in land vegetation cover, and changes in vegetation productivity over multiple years [6]. Satellite remote sensing data are widely employed for monitoring vegetation phenology across large geographical regions. In particular, the Moderate Resolution Imaging Spectroradiometer (MODIS) has emerged as a valuable tool due to its high temporal resolution (daily) and moderate spatial resolution (250–500 m). Vegetation indices, which consider various spectral signals, are valuable for quantifying vegetation dynamics [7,8]. The Ratio Vegetation Index (RVI), first proposed by Jordan in 1969 [9], provides relatively stable information for vegetation monitoring [10]. The Normalized Difference Vegetation Index (NDVI), developed through nonlinear processing of the RVI, stands as the most widely utilized vegetation index [11]. The NDVI effectively addresses irradiance biases associated with instrument calibration, solar zenith angle, terrain, cloud shadows, and atmospheric conditions, thereby enhancing its sensitivity to vegetation [12]. A significant improvement is the Enhanced Vegetation Index (EVI), introduced by Liu and Huete, which incorporates feedback terms to mitigate soil and atmospheric effects, thereby giving it an improved performance compared to the NDVI [13].

In vegetation phenology monitoring, the Normalized Difference Infrared/Water Index (NDII/NDWI) is commonly used to explain the impact of leaf moisture content [14]. Gonsamo et al. proposed the Phenology Index (PI), which combines the NDVI and NDII, to mitigate surface greenness effects, especially snow-mixing effects [15], significantly improving vegetation phenology monitoring [16]. Recently, Wang et al. introduced the Normalized Difference Phenology Index (NDPI), which incorporates RED, NIR, and SWIR bands, enhancing the spatial and temporal scalability of vegetation phenology monitoring [16]. Xu et al. conducted extensive field surveys in Inner Mongolia’s natural grasslands and demonstrated the NDPI’s superior spatial and temporal scalability for estimating aboveground fresh biomass (AGB) [17]. The NDPI is valuable for estimating the AGB in grasslands by reducing differences between soil backgrounds and distinguishing soil from vegetation. These advantages are particularly pronounced in grassland ecosystems, as the NDPI remains unaffected by snowmelt and exhibits sensitivity to initial vegetation growth during the greening period in spring, thereby enhancing the monitoring of spring phenology [18]. High-spectral-resolution data, such as hyperspectral data, are valuable for identifying herbaceous plant species within communities, enabling major species proportional area estimation and providing structural indicators for grassland ecosystems [19].

The distribution of grassland resources in Nanling is fragmented due to significant variations in elevation and steep slopes, resulting in low utilization rates [20]. The poor soil conditions lead to inferior grass quality compared to that of the northern grasslands [21]. Therefore, the key challenge in developing and utilizing grassland resources in mountainous areas is to protect the environment while ensuring sustainable grassland production [22]. This study aimed to fully capitalize on the advantages of the NDPI in grassland vegetation detection to conduct continuous phenological monitoring of the grassland ecosystem in the Nanling region and evaluate the change trend of the grassland biomass and the transition points of the time series in order to detect the spatial differentiation and time evolution of the southern mountain grassland ecosystem. To this end, we further explored the relationship between and thresholds of grassland coverage and NDPI values in the process of regional grassland landscape restoration, providing data support for ecological management oriented toward grassland restoration.

2. Materials and Methods

2.1. Study Area

The Nanling region spans parts of Hunan, Jiangxi, Guangdong, and Guangxi (Figure 1a). The area encompasses 4 provinces, 9 cities, and 70 counties, covering a total area of approximately 166,054 km² (Figure 1b). The predominant landforms in Nanling are hills and mountains, with an overall relative elevation difference of about 2000 m across the region [23]. Nanling serves as the watershed for two major river systems, namely the Pearl River and the Yangtze River, with a combined length of approximately 7750 km (Figure 1c). The Nanling mountain range significantly influences the local climate, acting as a barrier during winter to weaken the impact of dry and cold waves as they cross the ridge [24]. The Nanling region experiences an average annual precipitation of approximately 2037 mm and an average annual temperature of around 19.41 °C. The favorable hydrothermal conditions, extended grass period, and high productivity make Nanling a region with great production potential [25,26]. In order to facilitate the study and analysis, the study area in the Nanling region was divided based on the administrative divisions of the provinces and counties in China. Based on land use/land cover data, the identified grassland area in this region is about 11,000 km², accounting for approximately 6.6% of the total area (Figure 1d).

2.2. Methods

A comprehensive study was undertaken to examine the variations in NDPI values within the grasslands of the Nanling region (Figure 2). The investigation analyzed the interannual trends at 70 sites across the study area. Data were collected over a span of 22 years, from January 1st to December 27th each year, with the exception of missing data for the years 2000 and 2001. For each site, average NDPI values were calculated for roughly the same dates every year, creating a consistent time series of NDPI values spanning 22 years. To assess the trends in these time series, the Mann–Kendall (M-K) trend test was employed.

2.2.1. NDPI

The NDPI is defined by red (Band 1, 620~670 nm), NIR (Band 2, 841~876 nm), and SWIR (Band 6, 1628~1652 nm) bandwidths as follows:

$$\mathsf{\rho }=\mathsf{\alpha }\times \mathrm{RED}+(1-\mathsf{\alpha })\times \mathrm{SWIR}=0.74\times \mathrm{RED}+0.26\times \mathrm{SWIR}$$

(1)

$$\mathrm{N}\mathrm{D}\mathrm{P}\mathrm{I}={\displaystyle \frac{\mathrm{NIR}-\mathsf{\rho }}{\mathrm{NIR}+\mathsf{\rho }}}$$

(2)

where ρ is the weighted sum of the red and SWIR reflectance in order to minimize the difference between soil and snow, which are the main components of the land surface in the non-growing season, and α is a weighting coefficient that was set to 0.74 based on existing research [25]. In this study, the MODIS Surface Reflectance 8-Day 500 m product (MOD09A1), which is widely used in phenology monitoring studies [27,28,29,30], was utilized in TIMESAT 3.3 software to extract phenological parameters. The algorithm of the adaptive Savitzky–Golay filter employs a sliding window approach to linearly approximate the upper envelopes of vegetation index time series in a piecewise manner, gradually and accurately eliminating noise through iteration [30,31]. TIMESAT 3.3 includes the Seasonal Trend Decomposition by LOESS (STL) method, which utilizes the Locally Weighted Regression Smoother (LOESS) technique [32,33,34].

2.2.2. M-K Trend Test

The Mann–Kendall trend test method is commonly used in the study of vegetation dynamics and phenological responses [35,36,37,38,39,40]. The original hypothesis H₀ of the Mann–Kendall test is that in time series data (x₁, ..., x_n), there are n samples of independent and identically distributed random variables with no upward or downward trend. The alternative hypothesis H₁ is a bilateral test indicating a monotonic trend in the sequence, meaning that for any k and i (k ≠ i, i ≤ n), x_k and x_i have different distributions. The statistical variable S for the test is defined as follows:

$$\mathit{S}={\sum }_{\mathit{i}=1}^{\mathit{n}-1}{\sum }_{\mathit{k}=\mathit{i}+1}^{\mathit{n}}\mathit{sgn}$$

(3)

$$\mathit{sgn}({\mathit{x}}_{\mathit{k}}- {\mathit{x}}_{\mathit{i}})=\left\{\begin{array}{c}1,{ \mathit{x}}_{\mathit{k}}\text{>}{ \mathit{x}}_{\mathit{i}}\\ 0,{ \mathit{x}}_{\mathit{k} }={ \mathit{x}}_{\mathit{i}}\\ -1,{ \mathit{x}}_{\mathit{k}}\text{<}{ \mathit{x}}_{\mathit{i}}\end{array}\right.$$

(4)

where sgn denotes the sign function, and x_k and x_i are adjacent observations in the time series.

When n > 10, the statistical variable S is approximately normally distributed. If each number in the sequence is unique, the variance is

$$\mathrm{Var}(\mathit{S})={\displaystyle \frac{1}{18}}\mathit{n}\left(\mathit{n} - 1\right)\left(2\mathit{n}+5\right)$$

(5)

If there are non-unique data in the sequence, the variance is

$$\mathrm{Var}(\mathit{S})={\displaystyle \frac{1}{18}}\left[\mathit{n}\left(\mathit{n}- 1\right)\left(2\mathit{n}+5\right) - {\sum }_{\mathit{p}=1}^{\mathit{g}}{\mathit{t}}_{\mathit{p}}\left({\mathit{t}}_{\mathit{p}} - 1\right)\left(2{\mathit{t}}_{\mathit{p}}+5\right)\right]$$

(6)

In Equation (6), p is the number of repetitions, g is the number of unique numbers, and t_p is the number of repetitions for each repetition.

The statistic Z_MK for M-K trend test is

$${\mathit{Z}}_{\mathit{MK}}=\left\{\begin{array}{cc}{\displaystyle \frac{\mathit{S}-1}{\sqrt{\mathit{Var}\left(\mathit{S}\right)}}}, & \mathit{S}\text{>}0 \\ 0,& \mathit{S}=0\\ {\displaystyle \frac{\mathit{S}+1}{\sqrt{\mathit{Var}\left(\mathit{S}\right)}}},& \mathit{S}\text{<}0\end{array}\right.$$

(7)

For a given significance level α in a bilateral trend test, if |Z_MK| < Z_1−α/2, then the original assumption H₀ holds. If |Z_MK| ≥ Z_1−α/2, then the original assumption is unacceptable, indicating a significant monotonic trend in the sequence. If Z_MK > 0, it suggests an increasing trend in the sequence, and vice versa. When |Z_MK| ≥ 1.645, 1.960, 2.576, it indicates that the trend has passed the significance tests with confidence levels of 90%, 95%, and 99%, respectively.

In addition, the monotonic trend can be quantified using the Kendall slope β (Sen slope):

$$\mathsf{\beta }=\mathrm{median}\left({\displaystyle \frac{{\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{j}}}{\mathit{i}-\mathit{j}}}\right), 1 < j< i< n$$

(8)

In Equation (8), “medium” represents the median value. When β > 0, it indicates that the variable X has an upward trend; otherwise, it shows a downward trend [40].

The Mann–Kendall method can also detect the mutation points of time series. For time series x_i, construct a rank sequence r_ij to represent sample statistics for x_i > x_j (1 ≤ j ≤ i):

$${\mathit{r}}_{\mathit{ij}}=\left\{\begin{array}{cc}1,& {\mathit{x}}_{\mathit{i}}\text{>}{\mathit{x}}_{\mathit{j}}\\ 0,& {\mathit{x}}_{\mathit{i}}\text{≤}{\mathit{x}}_{\mathit{j}}\end{array}\right. 1\text{≤}\mathit{j}\text{≤}\mathit{i}$$

(9)

Define the sample cumulative S_k of order column r_ij as

$${\mathit{S}}_{\mathit{k}}={\sum }_{\mathit{i}=1}^{\mathit{k}}{\sum }_{\mathit{j}}^{\mathit{i}-1}{\mathit{r}}_{\mathit{ij}}, \mathit{k}=2,3,\text{…},\mathit{n}$$

(10)

Assuming that the time series is random and independent, UF_k, the statistic of the standard normal distribution, is defined as

$${\mathit{UF}}_{\mathit{k}}={\displaystyle \frac{{\mathit{S}}_{\mathit{k}}-\mathit{E}\left({\mathit{S}}_{\mathit{k}}\right)}{\sqrt{\mathit{Var}\left({\mathit{S}}_{\mathit{k}}\right)}}}, \mathit{k}=2,3,\text{…},\mathit{n}$$

(11)

In Equation (11), E (S_k) is the mean of S_k and Var (S_k) is the variance of S_k.

$$\mathit{E}\left({\mathit{S}}_{\mathit{k}}\right)={\displaystyle \frac{\mathit{k}(\mathit{k}+1)}{4}}$$

(12)

$$\mathit{Var}({\mathit{S}}_{\mathit{k}})={\displaystyle \frac{\mathit{k}(\mathit{k}-1\left)\right(2\mathit{k}+5)}{72}}$$

(13)

At a given significance level α, if |UF_k| > U_α/2, it indicates a significant trend change in the sequence. Using the inverse UB_k of the time series x_i, repeat the above process. Through further analysis of the statistical sequences UF_k and UB_k, the time nodes of sequence x_i mutations can be obtained, and the mutation regions of the sequence can be displayed. If UF_k > 0, it suggests an increasing trend in the sequence; otherwise, it shows a decreasing trend. When UF_k and UB_k exceed any critical line, it signifies a significant trend of an increase or decrease in the sequence. If the curves of the UF_k and UB_k sequences intersect between critical lines, the moment of intersection marks the moment of mutation.

2.2.3. Segmented Regression

After the mutation point (date) of the NDPI time series was identified by Mann–Kendall method, the maximum, average, and minimum NDPI values corresponding to the mutation date were obtained. We then conducted segmented linear regression between the grassland coverage of each county and the values of the maximum, average, and minimum NDPI values at the change points. We used the “segmented” package in R 4.3.2 to fit the segmented linear regression models described above, and then detected breakpoints for these models [41].

2.3. Data Sources

The land use/land cover data were obtained from the Resources and Environmental Sciences and Data Center (https://www.resdc.cn/, accessed on 19 March 2023), covering the years 2000, 2005, 2010, 2015, and 2020, with a spatial resolution of 30 m. Patches classified as grassland in the land use/land cover data were used as grassland boundaries for this study. The grass cover areas that persisted for at least 3 out of the 5 time points (2000 to 2020) were identified as relatively stable grassland areas in the Nanling region, and other ubefore 25 January 2025 nstable grasslands were excluded. A total of 70 counties (districts) were taken as statistical samples, and for each county (district), the mean, maximum, and minimum NDPI values of grassland patches were computed.

The MODIS surface reflectance 8-Day 500 m images (MOD09A1) were obtained from the National Aeronautics and Space Administration (https://ladsweb.modaps.eosdis.nasa.gov, accessed on 19 March 2023). A total of 1004 MODIS images of the Nanling region were acquired, covering the period from 2000 to 2021. The continuous time series of the NDPI values of the grasslands in the Nanling region was calculated in Google Earth Engine (GEE). The 8-day product can effectively reduce cloud contamination [17], but to avoid cloud interference, we only kept those pixels without clouds, cloud shadows, and cirrus clouds. In TIMESAT 3.3, the original sequence of the NDPI grid means was reconstructed using the Savitzky–Golay filter with the sliding window width set to 5.

3. Results

3.1. Grass Recovery Pattern During 2000–2021

There was a significant trend in the recovery of grass in the Nanling region, and 2005 marked the turning point of grassland recovery during the past two decades (2000–2021) (Figure 3). The average NDPI values exhibited a significant increasing trend. This trend was more pronounced in the southeast compared to the northwest (Figure 4a), aligning with the distribution of annual average temperature in the Nanling region. The increase in NDPI values was most noticeable around 2010 to 2014, with a sequential pattern observed in the eastern, western, and central regions (Figure 4b).

The average Z_MK statistic for the M-K test in each county (or district) within the Nanling region was approximately 9.198, which exceeds the significance threshold of 99.9% (Z_MK > 3.239), indicating a significant increasing trend in the NDPI values across different sites from 2000 to 2021 (Figure 4a). Overall, there was a clear and consistent pattern of increasing grassland biomass in the Nanling region, implying an improvement in grassland quality and the absence of overall degradation over 2000–2021. However, there is also a spatial difference in the significance of the NDPI increasing trend, and the significance level in the southeast of Nanling was higher than that in the northwest (Figure 4a), which highlights the substantial influence of thermal conditions on grassland productivity within the region. The Z_MK values exceeding 11 were predominantly concentrated in Guangdong (mean Z_MK = 10.804) and Jiangxi (mean Z_MK = 10.627). In Guangxi (mean Z_MK = 8.148) and Hunan (mean Z_MK = 7.756), the significance of the NDPI growth trend was relatively low. Among all the sample sites, the NDPI trend was least significant in northern Guangxi, especially in Ziyuan County (Z_MK = 4.242).

The M-K trend test revealed relatively small magnitudes for the slope parameter β for the NDPI values due to the long time series. The range of β values across the different sample sites in the Nanling region was 4.8×10⁻⁵ to 11.0 × 10⁻⁵, with an average value of approximately 7.7 × 10⁻⁵. The eastern and southern parts of the region exhibited larger growth rates in the NDPI values, while the western region showed smaller rates (Figure 4a). The average β for the NDPI changes across the different provinces/regions was as follows: Jiangxi (8.02 × 10⁻⁵) > Guangdong (7.97 × 10⁻⁵) > Hunan (7.58 × 10⁻⁵) > Guangxi (7.26 × 10⁻⁵). The statistical parameter Z_MK in the M-K trend test showed good consistency with the slope parameter β. In the regions where Z_MK had a higher significance, the corresponding NDPI increasing trend tended to have relatively larger slope values. However, the spatial differentiation pattern of the overall NDPI trend (β) in the Nanling region was not as distinct as the regional patterns observed regarding the significance of Z_MK.

The sequences of the mean NDPI value of the grassland in each county experienced a turning point mainly between 2010 and 2014 (Figure 4b), indicating that these five years are the main period for grassland recovery from disturbance. Among all the counties, Ziyuan, Chengbu, and Jiahe stood out with the earliest change points in September and October 2000, which significantly deviated from the other counties. These sites also exhibited the lowest significance in their NDPI increasing trends (Z_MK < 5.1), suggesting relatively stable NDPI variations and a less pronounced increase in grassland biomass over the study period. In Jiangxi province, the majority of the sites experienced change points in 2010 and 2011, while in the Guangxi region, it mainly occurred in 2012 and 2013. In the central region, the Guangdong and Hunan provinces began to display an upward trend between 2013 and 2015. The latest change points were situated in the central-eastern part of the Nanling region.

3.2. Grassland Recovery Monthly Threshold

The annual variation in NDPI values showed an increasing trend, indicating improved productivity and reduced degradation during different seasons. The months of January to March and June to September were identified as the threshold months for these changes. The peak of the average NDPI value primarily occurred from June to September (Figure 5a), aligning with the peak growing season of grasslands. Throughout the 22-year study period, June, July, August, and September were the months with the highest frequency of peak values, occurring 8, 5, 7, and 2 times, respectively. For example, in 2001, 2003, and 2016, nearly 50% of the regions peaked in July, more than 70% in June, and more than 50% in August, respectively. In some areas, the peak values occurred in April or May, which may suggest unfavorable weather conditions such as flooding or hot and dry conditions during the typical peak season of June to September. These conditions could impede grassland biomass accumulation, resulting in an earlier peak in the NDPI values. Conversely, peaks occurring in October may be attributed to an extended growing season with favorable water and thermal conditions, allowing continued grassland growth and organic matter accumulation.

The minimum NDPI values were primarily observed in January to March and November to December (Figure 5b). Among the 22 years, February had the highest frequency of being the month with the minimum NDPI value in 15 years, while January and March had the highest frequency in the remaining years. Most grasslands reached their minimum aboveground biomass at the end of winter, resulting in the lowest NDPI values primarily occurring in February. The NDPI value reaching its minimum in March may indicate prolonged periods of dry and cold climates. Conversely, the minimum NDPI values in January may be influenced by early warming and an advancement in the grassland growing season.

The upper limit of the average grassland NDPI value in the Nanling region ranged from 0.572 to 0.727 (Figure 6a). The highest upper limit of the NDPI variation were found in Heping and Liannan in Guangdong, Dinnan in Jiangxi, and Longsheng in Guangxi, indicating that these grasslands have greater potential productivity. Conversely, the lowest upper limits were observed in Ganxian and Yudu in Jiangxi, Jiahe in Hunan, and Longchuan in Guangdong, where the NDPI upper limits fell below 0.6 (Figure 6b), mainly in areas with relatively low precipitation and a dense population distribution.

3.3. Seasonal Characteristics of Grassland NDPI Changes

Using the STL decomposition method, the length of the growing season and its start and end dates were extracted for the various grassland sites (Figure 7a). The length of the growing season varied from 188 to 247 days, with longer growing seasons observed in the southern regions compared to the northern regions. The interannual variation in the start and end dates of the growing season was minimal, with the deviations within a span of 5 days. The range of start dates for the growing season spanned from approximately March 29th to May 9th, with the average start date ranging between April 7th and May 1st. Notably, around half of the sites commenced their growing season on April 15th, while the site with the latest start date was located in the northeastern part of the study area in Jiangxi. In terms of the end date of the growing season, it typically fell between October 24th and December 27th. However, due to the influence of mountain ranges and variations in locations on the north and south slopes, significant differences were observed between some sites. Generally, there was a trend in the later end dates in the southern region compared to the northern region.

The grasslands in Guangdong and Guangxi displayed early start dates and late end dates for the growing season, lasting for more than 210 days. This pattern was strongly correlated with the higher precipitation and relatively higher temperatures experienced in these regions. Conversely, the grassland sites with growing seasons lasting less than 200 days were predominantly located in Jiangxi and Hunan provinces, specifically on the northern slope of the Nanling Mountains. In the northeastern region, there was a delay in the start of the growing season and an earlier end, which could be attributed to the delayed spring warming and early cooling in autumn and winter in these areas. The shorter growing season also set an upper threshold for the NDPI variation, with some sites exhibiting lower NDPI values during the winter degradation.

We derived the average amplitude of the annual variation in the mean NDPI values (Figure 7b). The results reveal a distinct northwest-to-southeast gradient, which exhibited an inverse correlation with the spatial distribution of the annual average temperatures in the Nanling region. Specifically, sites with higher temperatures tended to display smaller NDPI amplitudes, suggesting that insufficient heat may contribute to significant fluctuations in grassland productivity. In contrast, areas with ample heat experienced more stable growth and maintained higher biomass levels, resulting in relatively weaker degradation during the winter season (Figure 7b). The eastern regions of Jiangxi and Guangdong encompassed the most stable grasslands, while Hunan exhibited the largest amplitude of variation. This pattern aligned with the lower precipitation levels in Hunan compared to the Guangxi region, indicating that the stability of grasslands in the western areas was additionally influenced by water availability. Consequently, the average amplitude of the NDPI variation in the Nanling region was primarily controlled by temperature, with the western areas further affected by precipitation constraints.

3.4. Peak Value Changes in Grass NDPI

The maximum NDPI values indicated an improvement in grassland productivity in the Nanling region. It is crucial to delve deeper into the statistical analysis of the annual NDPI values to gain a better understanding of the potential for grassland biomass recovery and the limits of degradation. By examining the annual NDPI values, we can assess the variations and fluctuations in the grassland biomass. This information will help identify the resilience of grasslands and their ability to recover from degradation, as well as provide insights into the factors that contribute to the limits of degradation.

From 2000 to 2021, the highest NDPI values across all the sites were primarily observed from May to September, accounting for 96%. Among these months, the period from June to August represented approximately 79% (Figure 8a). The maximum NDPI values of different sites in the Nanling region exhibited an increasing trend, with the majority ranging from 0.7 to 0.8, indicating an improvement in grassland productivity over the past 20 years (Figure 8b). In particular, the maximum NDPI values in Guangdong and Guangxi provinces were relatively higher compared to other areas. This can be attributed to the abundant precipitation brought by the maritime moisture received on the southern slope of the Nanling Mountains. The presence of this moisture creates favorable conditions for grassland growth, allowing for the accumulation of aboveground biomass. However, in certain areas, such as Yizhang, Yongxing, Jiahe, Linwu, and Xintian in Hunan Province, the growth potential of the NDPI values was significantly lower. These five counties are situated on the northern side of Mengzuling, Qitianling, and other surrounding mountains. This limitation in moisture availability hampers the development of grassland productivity in these specific regions.

From 2000 to 2021, the minimum NDPI values within the Nanling region primarily occurred in January, February, and December, accounting for approximately 67.7% of the minimum NDPI values (Figure 8c). Among these months, January represented approximately 19.3%, while February represented approximately 24.4%. The minimum NDPI values were observed to be between 0 and 0.14 (Figure 8d). Several sample sites exhibited minimum NDPI grid values below 0.02, indicating a relatively low level of vegetation cover approaching completely bare soil. Examples of such sites include Ganzhou, Yudu, Yongzhou, Ruyuan, Lianzhou, Quanzhou, and Fuchuan. These sites are predominantly located in hilly plains or mountain passes between the mountains on the northern slopes of the Nanling. Consequently, sample sites situated on the southern side of the mountains, such as Anyuan, Dingnan, Quannan, Lianshan, and Liannan, experienced relatively less severe local degradation of grassland. Overall, the variation in the minimum NDPI values among different sites within the Nanling region can be attributed to a combination of factors, including local weather conditions, topography, and human activities.

3.5. NDPI Thresholds for Mountain Grassland Restoration

In 63 out of 70 counties, the maximum, mean, and minimum NDPI values corresponding to the time series turning points could be detected. We constructed piecewise linear regression models for the grassland coverage and the maximum, mean, and minimum NDPI values corresponding to the time series turning points in 63 counties (Figure 9). The results showed that the mean NDPI values were significantly (p < 0.05) correlated with grassland coverage (Figure 9a), but the maximum and minimum NDPI values were not significantly (p > 0.05) correlated with grassland coverage (Figure 9b,c). When the grassland coverage ranged from 0 to 2.1, the mean NDPI value of the grassland patches corresponding to the grassland productivity recovery gradually increased to 0.5; however, when the grassland coverage was greater than 2.1, the mean NDPI value of the grassland patches corresponding to the grassland productivity recovery still fluctuated around 0.5 (Figure 9a). This means that once the grassland coverage rate in the Nanling region exceeds 2.1%, the mean NDPI value corresponding to grassland recovery will not vary significantly. This indicates that as long as the grassland coverage in the Nanling region is maintained above 2.1%, the restoration efficiency of the grassland landscape can be considered satisfactory. The grass landscape restoration in the counties of the Nanling region might be mainly characterized by a low level of disturbance from human activity during the past two decades.

4. Discussion

4.1. Ecological Restoration of Grassland in Southern Mountainous Area of China

Remote sensing monitoring of grasslands has been used in recent years to assess large-scale changes in grass landscapes. Many remote sensing-related indicators have been identified as effective in detecting shifts in grassland degradation and restoration [42]. Previous studies have extensively examined the impact of soil conditions and moisture sensitivity on remote sensing monitoring of grasslands [43,44,45]. These studies have demonstrated the considerable advantage of the NDPI in estimating aboveground biomass through model evaluations [17]. However, it is important to note that there are substantial differences between southern grasslands and northern prairie ecosystems in China. The accuracy of the NDPI in biomass estimation within the Nanling grasslands, as well as its superiority over other indices such as the NDVI and EVI in biomass models [46,47,48], lacks validation from ground survey measurements. Therefore, further comparisons and evaluations are necessary to evaluate the consistency and applicability of the NDPI in monitoring changes in southern mountainous grassland ecosystems. Additionally, it would be valuable to establish a relationship between biomass and productivity by considering different grassland types and hydrothermal conditions [49]. This could involve the development of a remote sensing estimation model for the Net Primary Productivity (NPP) in southern grasslands based on the NDPI [50], thereby broadening the applications of the NDPI in this region [51]. To validate the suitability and effectiveness of the NDPI in monitoring changes in southern mountainous grassland ecosystems, we recommend conducting comparative studies and assessments that incorporate ground-truth data and field investigations [52]. These efforts will enhance our understanding of the consistency and applicability of the NDPI in monitoring the dynamics of southern grassland ecosystems [53].

In this study, the investigation into the influencing factors of the NDPI was limited, and the mechanisms behind various potential factors remain unclear. The dominant roles of precipitation and temperature in the different regions, as well as the extent of the impact of human activities, have yet to be determined. Specifically, further research is needed to explore the correlation between the NDPI trend changes over years and solar activity cycles. It is important to comprehensively consider factors such as elevation, slope, rainfall, temperature [54,55,56], land use/land cover changes, grassland types, and human disturbances to examine their combined influence on grassland degradation and productivity development [57,58]. We fitted the linear relationship between the multi-year mean NDPI value series of the grassland in each county and the spatial mean value series of the four impact factors (Figure 10). The results show that the NDPI value was significantly (p < 0.001) negatively correlated with population density (Figure 10a), positively correlated with precipitation and terrain slope (Figure 10b,d), and had no significant linear relationship with temperature (Figure 10c). This indicates that the grassland phenology in Nanling area is controlled by both human activities and natural environmental conditions. Human activities may mainly influence grassland coverage changes under varying levels of human disturbance intensity, while natural environmental conditions may primarily affect grassland coverage regreening.

4.2. Limitations and Perspective

There are several limitations related to the data and methods employed in this study. Specifically, the use of the GlobeLand30 dataset as a substitute for the unavailable 2020 land use data introduces inconsistency in the data sources, which could potentially affect the results [59,60]. Although efforts were made to minimize the impact by extracting samples from grids identified as grassland across multiple periods, striving for consistency in data sources for future studies would enhance the robustness and reliability of the findings [61,62,63].

Furthermore, the fragmented nature of the grassland landscape in the Nanling mountainous area poses a challenge for accurately capturing the dynamics of grassland growth using 500 m resolution MODIS imagery [64,65]. To improve spatial resolution, exploring alternative approaches like spatial–temporal fusion to resample the data to a finer resolution, such as 30 m, would enable a more detailed investigation of the spatiotemporal distribution characteristics of the NDPI values. Additionally, considering the extraction of similar bands from Landsat imagery for comparative analysis could help assess the influence of different image sources on the performance and superiority of the NDPI [66].

When conducting data smoothing reconstruction and phenology extraction using TIMESAT, it is important to acknowledge that the selection of experimental methods and parameter settings is subjective [67]. Insufficient testing and comparisons of various processing approaches were conducted, which could potentially affect the accuracy of the reconstruction and extraction results [68]. It should be emphasized that the extracted growing season should be considered as a reference rather than an absolute measure. To overcome this limitation and strengthen the analysis, we recommend conducting comprehensive sensitivity analyses that explore a range of methodological choices and parameter settings. This approach will provide a deeper understanding of the uncertainties and variations associated with different approaches, thereby facilitating the selection of the most suitable methods for data reconstruction and phenology extraction [69]. Moreover, the inclusion of ground truth data or independent validation datasets is crucial for assessing the accuracy and reliability of the results.

Overall, this analysis offers critical insights into the trends, thresholds, and spatial patterns of grassland ecosystem restoration in the Nanling region. The results indicate the net effects of climatic factors and human activities on grassland productivity. These findings contribute to a better understanding of the dynamics and characteristics of grassland ecosystems in the region, which can provide guidance for agricultural development, land management, and grassland ecosystem conservation in the Nanling region. This knowledge is essential for promoting the sustainable development of grassland ecosystems and advancing the integrated protection of ecological security barriers.

5. Conclusions

The grasslands in the hilly areas of the middle and low mountains in southern China have high restoration efficiency and ecological value due to the favorable combination of water and heat conditions [70]. However, our understanding of the dynamic adaptation process of grassland restoration under the combine effects of climate change and human activities remains limited. Utilizing remote sensing technology to monitor its phenological trends can help understand the adaptive characteristics of grassland restoration, and the restoration threshold is a key turning point for identifying stages of change in grassland restoration. Therefore, based on the widely distributed grasslands in the region, we employed the NDPI to monitor their recovery process and applied statistical regression methods to detect their change thresholds. This study yielded insights into overall trends, seasonal variations, and their thresholds and spatial characteristics. The conclusions were as follows. (1) Overall increasing trend: The average NDPI values exhibited a significant increasing trend. This trend was more pronounced in the southeast compared to the northwest, aligning with the distribution of annual average temperature in the Nanling region. The most notable increase in the NDPI values occurred between 2010 and 2014, with a sequential pattern observed in the eastern, western, and central regions. (2) Seasonal characteristics and thresholds: The annual variation range of the NDPI values showed an increasing trend, suggesting improved productivity and reduced degradation during different seasons. January to March and June to September were identified as the threshold months for these changes. (3) Growing season length: The length of the growing season varied from 188 to 247 days, with longer growing seasons observed in the southern regions compared to the northern regions. (4). Spatial differentiation in productivity and degradation: The analysis of the NDPI grid maximum values indicated an improvement in grassland productivity in the Nanling region, with better growth conditions observed in Guangdong and Guangxi provinces. (5) Mountain grass landscape restoration was closely related to grassland coverage and the mean NDPI values in grassland patches during the past two decades. When the grassland coverage was greater than 2.1%, the grassland landscape restoration in the region was satisfactory.

The analysis of the grassland NDPI values in the Nanling region from 2000 to 2021 revealed an increasing trend and an overall improvement in productivity potential. The spatial distribution of grassland productivity follows a pattern of higher values in the southeast and lower values in the northwest. This distribution is influenced by temperature, precipitation, and human activities. To ensure sustainable development and conservation of grassland ecosystems, targeted management strategies should be implemented, particularly in regions where human factors contribute to significant fluctuations in grassland productivity. These strategies should consider local climate characteristics and their impact on grassland biomass. By leveraging favorable conditions such as water availability and heat, and addressing challenges such as topography and grass quality, efforts can be focused on the development of grassland livestock production and the implementation of supportive policies. Increased investment in technology and the establishment of efficient industrial development models are crucial for the growth of grassland livestock in the southern region.