Vegetation dynamic analysis based on multisource remote sensing data in the east margin of the Qinghai-Tibet Plateau, China

Image sources

Data processing

(1) Image sampling.

Sampling is indispensable when finding the trend analysis of NDVI at the site scale. The location and the size of the sample especially will impact the results of NDVI trend analysis. Many factors were considered, such as the sample being an area of vegetation cover, and only using the single type of vegetation. Moreover, mixed pixels were eliminated for NDVI data with two different resolutions. Using the analysis above, we selected 15 representative sample sites, with a sample size of 8 km × 8 km (GIMMS 3g NDVI, 1 × 1 pixel. VGT\PROBA-V NDVI, 8 × 8 pixels). The 15 sample sites selected were Songpan (SP), Kangding (KD), Litang (LT), Jiuzhi (JZ), Derong (DR), Nangqian (NQ), Tongren (TR), Dege (Dege), Ruoergai (REG), Maerkang (MEK), Ruobei (RB), Zeku (ZK), Gongshan (GS), and Yanyuan (YY). The location and vegetation types of each site are shown in Table 2 and Fig. 2.

Table 2:

Vegetation types and their use of the sample site.

Site	Longitude (°E)	Latitude (°N)	Altitude (m)	Vegetation type	Data use
Songpan (SP)	103.567	32.650	3732	Needle-leave deciduous forest (NEF)	Phenology analysis and vegetation dynamic
Ruobei (RB)	103.894	34.105	3155	Broadleaved deciduous forest (BDF)	Vegetation dynamic
Zeku (ZK)	101.504	35.094	3889
Gongshan (GS)	98.667	27.750	2297
Yanyuan (YY)	101.517	27.433	2571
Kangding (KD)	101.967	30.050	4038	Broadleaved evergreen forest (BEF)	Vegetation dynamic
Litang (LT)	100.267	30.000	4581
Jiuzhi (JZ)	101.483	33.433	3985	Alpine meadow (AM)	Phenology analysis and vegetation dynamic
Drerong (DR)	99.650	33.750	4396
Nangqian (NQ)	96.483	32.20	3661
Tongren (TR)	103.389	33.796	3164	Typical meadow (TM)	Vegetation dynamic
Dege (DG)	98.583	31.800	3918
Ruoergai (REG)	102.967	33.583	3444	Plain grassland (PG)	Vegetation dynamic
Maerkang (MEK)	102.451	31.896	4708	Swamp (SW)	Vegetation dynamic
Ruoe (RE)	102.017	35.517	3229	Farm land (FL)	Phenology analysis and vegetation dynamic

DOI: 10.7717/peerj.8223/table-2

(2) Sample statistics.

The basic principle of the sample statistics (Fig. 3) is that the calculation of pixel parameters was performed by setting a zonal window with a fixed size. The calculated parameters included the sum, mean, maximum, and minimum of the pixel values in the zonal window. The size of the zonal window was set to 8 km × 8 km according to the lower resolution of the two different NDVI images in this study. The results were divided into three cases (possibilities) for GIMMS 3g images, using the mean statistics of a GIMMS 3g image as an example. The formulations are as follows.

(2)${\displaystyle {\text{Case}}_1:\text{Mean}=V_{1}or{V}_{2}or{V}_{3}or{V}_4,}$

(3)${\displaystyle {\text{Case}}_2:\text{Mean}=\frac{V_1+V_2}{2}or\frac{V_1+V_3}{2}or\frac{V_2+V_4}{2}or\frac{V_3+V_4}{2},}$

(4)${\displaystyle {\text{Case}}_3:\text{Mean}=\frac{V_1+V_2+V_3+V_4}{4},}$

Where V is the pixel value of GIMMS 3g NDVI images. According to this principle, the formula above was additionally followed for VGT/PROBA-V NDVI image statistics. Due to its higher resolution, more pixels were counted than in GIMMS NDVI images. Compared to the Point to Value Method (PVM), this method provides more accurate statistical results and reduces invalid values.

Detection of phenological change

Previous studies (Chandola et al., 2010; Chu et al., 2018; Gu et al., 2018; Tang et al., 2018) have used the dynamic threshold method to define the start of the growing season (SGS), the end of the growing season (EGS), and the length of the growing season (LGS). This method defines SGS and EGS as the moment when the pixel values increase and decrease, respectively, to a certain proportion of NDVI amplitude in one year. This method can avoid the mutual interference caused by different regional water-heat conditions, vegetation, and soil types in the fixed threshold method. The specific formulation is as follows:

(5)${\displaystyle D_{SGS,EGS}=\left(NDV{I}_{max}-NDV{I}_{min}\right)\times C,}$

Where NDVI_max is the maximum value and NDVI_min is the minimum value of NDVI in one year. C is the threshold ratio. According to previous studies (Yu, Luedeling & Xu, 2010; Jonsson & Eklundh, 2002) and the annual variation of vegetation in the study area, the SGS and EGS thresholds were set at 20% in the sample site. D is a time judgment function where the time point is taken as EGS or SGS when NDVI increases or decreases to the threshold C. The time span from SGS to EGS is the LGS.

Trend analysis of the NDVI change

The stability and accuracy of the trend analysis method have been tested and it is widely used (Chu et al., 2018; Gu et al., 2018; Tang et al., 2018). This method was applied at a pixel scale in the period 1982–2018 to detect NDVI change and distribution. The formula is:

(6)${\displaystyle slope=\frac{n\times \sum _{i=1}^n\left(i\times NDV{I}_i\right)-\sum _{i=1}^{n}i\sum _{i=1}^{n}NDV{I}_i}{n\times \sum _{i=1}^n{i}^2-{\left(\sum _{i=1}^{n}i\right)}^2},}$

where i is the number of the year, with the range 1 to n (GIMMS NDVI, n = 34; VGT NDVI, n = 20). NDVI_i represents the growing season NDVI of year i. The slope is the trend in NDVI from 1982 to 2018. A slope >0 indicates an increasing trend, while a slope <0 represents a decreasing trend during the study period.

Mann-Kendall trend test

The Mann-Kendall (M-K) trend test has been widely used as a non-parametric statistical test method (Gocic & Trajkovic, 2013; Yavuz, 2018). This test method does not require samples for normal distribution, and it has a high degree of quantification. This method is suitable to trend test time series data, especially the data of long-term series. The principles of the M-K trend test are as follows.

(7)${\displaystyle S_k=\sum _{i=1}^k{r}_i,k=2,3\cdots ,n}$ (8)${\displaystyle r_i=\left\{\begin{array}{c}+1x_i>x_j\hfill \\ +0x_i<x_j\hfill \end{array}\right.j=1,2\cdots ,i}$

Where n is the number of the sample, and x is the time span. s_k is the cumulative number of values at the time i,  greater than those at the time j.

(9)${\displaystyle U{F}_k=\frac{\left|S_k-E\left(S_k\right)\right|}{\sqrt{var\left(S_k\right)}},k=1,2\cdots ,n}$ (10)${\displaystyle \left\{\begin{array}{c}E\left(S_k\right)=\frac{k\left(k-1\right)}{4}\hfill \\ var\left(S_k\right)=\frac{k\left(k-1\right)\left(2k+5\right)}{72}\hfill \end{array}\right.,k=2,3,\dots ,n}$

Where UF_k is a statistic calculated by x₁, x₂, …, x_n. $E\left(S_k\right),$ and var(S_k) is the mean and variance of S_k.UB_k= −UF_k(k = n, n − 1, …, 1), UB₁=0. a is the signification level. If a = 0.05,  then the u_0.05 =  ± 1.96.

A UF_k > 0 indicates an increasing trend, while a UF_k < 0 represents a decreasing trend. When the critical line is exceeded (u_0.05 =  ± 1.96.), a significant upward or downward trend is shown.

Bayesian trend test

The traditional approach of trend analysis cannot detect the turning point for a set of data in a time series, and cannot estimate the change after it happens. The Bayesian method can overcome these difficulties. The Bayesian theorem constructed the probability distribution (posterior probability) function of the turning point of data change trends. The specific formula is:

(11)${\displaystyle P\left(\theta |\text{data}\right)=\frac{P\left(\text{data}|\theta \right)\times P\left(\theta \right)}{P\left(\text{data}\right)},}$

Where P(θ) denotes the prior probability distribution. P(θjdata) is the posterior probability distribution, representing the confidence level of a turning point of change trends. P(data) indicates the evidence; more importantly, it indicates the moment of data peak. For example, in a set of data in time series X = {x₁₉₈₂, x₁₉₈₃, x₁₉₈₄ ⋯ x₂₀₁₅}, how do the data change? We would need to calculate the moment and the probability of the turning point occurring when the trend changes. We referred to previous studies (Hèou et al., 2017; Peter et al., 2018; Zhao et al., 2019; John, Thomas & Christopher, 2015) to calculate the moment and probability of the turning point occurring (Fig. 4). The combination of the Bayesian approach and M-K trend analysis can improve the accuracy of trend analysis, and can validate M-K trend analysis. The flowchart of methodology used in this study is shown in Fig. 5.

Phenological trends of vegetation

The dynamic threshold method was used to detect the phenological changes of three ecosystems from 1982 to 2015. The trend of phenology change (Figs. 6A, 6B, 6C, 6G, 6H, and 6I) and the M-K trend test (Figs. 6D, 6E, 6F, 6J, 6K, and 6L) are shown in Fig. 6. In the forest ecosystem (6A and 6E), the SGS change trend rose (slope = 10.08 days/decade, UF <0, Sig. = 0.05), and then began to fall (slope = 2.06 days/decade, UF >0, Sig. = 0.05), showing overall trending delays (slope = 0.66 days/decade). The EGS showed a slight decline from 1982 to 2015 (slope = 0.15 day/decade), although it began to advance after 1998 (slope = 0.89 day/decade). It was calculated that the LGS decreased by 0.51 days/decade based on the SGS and EGS during the study period. The curves of the SGS and EGS for grassland ecosystems are shown in Figs. 6C and 6I. The SGS of the grassland ecosystem was accelerated by 0.64 days/decade, and the EGS was delayed by 0.46 days/ decade. Accordingly, the LGS was prolonged by 1.1 days/decade from 1982 to 2015. The variation trend of SGS in the grassland ecosystem was consistent with that of the forest ecosystem. The year 1998 was the turning point of the curve, and after that year, the SGS curve showed trending delays (slope = 2.68 days/decade). For the farmland ecosystem (Figs. 6A and 6G), the trend of change was consistent with that of the grassland and forest ecosystem, except for the change rate of SGS. From 1982 to 2015, the SGS trend was accelerated by 2.98 days/decade, and the EGS was delayed by 1.07 days/decade. Due to the acceleration of the SGS and the delay of the EGS, the LGS was prolonged by 8.74 days/decade (1982-2015).

Spatial dynamics of vegetation

We calculated the trend and spatial distribution of NDVI changes across different seasons from 1982 to 2015 (Fig. 7). In spring (Fig. 7A), the trend range of NDVI change was from -0.0679/decade to 0.0693/decade. Vegetation increased most in the area with a percentage of 82.1%, and were mainly distributed in SP, GS, YY, REG and MEK. The needle-leaf deciduous forest showed vigorous activity in SP, GS, and YY, and the vegetation coverage of swamp and plain grassland also experienced an upward trend in REG and MEK. The vegetation dynamics showed a weak upward trend in DR, JZ, and NQ. The slope was 0.005/decade to 0.013/decade in the alpine and subalpine meadows distributed in this area. In addition, vegetation deterioration had occurred in regions of the EMQTP, such as the west of ZK and DR. The west of DR was especially affected, as it is a desert grassland with low vegetation coverage. As shown in Fig. 7B, summer vegetation experienced great degradation with more than 51.5% of the area showing a downward trend. In particular, the south section of the study area had a slope of 0.08 / decade. The proportion of the vegetation increase area was 71.7%, and was mainly distributed in NQ and DG regions due to the degeneration of desert grassland. The vegetation dynamic in autumn (Fig. 7C) was basically consistent with the mean trend of the full growing season (Fig. 7D). The area of NDVI increase accounted for 75.1% in the full growing season. The difference between autumn and the growing season is only shown in the slope rate of NDVI increase and decrease. In the growing season, the slope was -0.0394/decade to 0.0486/decade. The maximum value appeared in RE, the south of SP, and the central region of the study area. Through the analysis of vegetation dynamics across different seasons, we found that the decrease of NDVI was mainly distributed in the western DR region and near the dry valley in the south of the study area. The areas with obvious increases were distributed in the east section of the study area. Because the EMQTP is a sub area of QTP, the spatial distribution of vegetation change provided a quantitative value that could be compared with Peng’s results (Fig. 7 of Peng’s work). Although the time ranges of NDVI were greatly different (1982-2003 and 1982-2015), they can better reflect the decreasing trend of vegetation in the arid valley.

Phenological trends of vegetation

We used the VGT\PROBA-V NDVI with a fine resolution to calculate phenological trends. The change detection of the SGS and EGS were performed on the most recent 19 years. The trend of SGS change for the forest ecosystem showed delay (slope = 0.66 days/decade, −1.96 <Sig. <1.96), and the EGS showed a slight advanced trend (slope = 2.99 day/decade, UF <0, −1.96 <Sig. <1.96). Thus, it was calculated that the LGS decreased by 3.65 days/decade during this period. The SGS and EGS curves for the grassland ecosystem are shown in Figs. 8C and 8G. The SGS of the grassland ecosystem was delayed by 0.92 days/decade, and the EGS was delayed by 1.05 days/ decade. Accordingly, the LGS was prolonged by 0.03 days/decade. Compared with the forest ecosystem, the grassland ecosystem showed an obvious uptrend. In the farmland ecosystem (Figs. 8A and 8D), the SGS was delayed by 5.76 days/decade and the EGS was delayed by 0.67 days/decade. Due to the delays of the SGS and the EGS, the LGS was shortened by 5.09 days/decade.

Spatial dynamics of vegetation

In the most recent 19 years (1999-2018), NDVI showed a decreasing trend with the area accounting for 63.4% in spring (Fig. 9A). Most was grassland cover area with a slope less than 0.004/decade. Moreover, NDVI had decreased significantly in the eastern and southern dry valleys of the study area. This area was covered by forest, especially obvious in the east of REG, south MEK, and the northwestern section of KD. The areas of NDVI increase in summer (Fig. 9B), autumn (Fig. 9C), and the growing season (Fig. 9D) were 71.8%, 60.9% and 65.4%, respectively. The slopes of NDVI increase and decrease were similar and the distribution area was roughly consistent. Compared with the spatial distribution of NDVI changes from 1982 to 2015, a large area near REG showed some obvious differences. This area was covered by alpine and sub-alpine meadows, and NDVI showed a downward trend with a slope of about 0.005/decade, especially in the summer (Fig. 9B). This indicated that vegetation degradation in the NEG region has been more noticeable in the last 19 years, and more attention should be paid to the changes in the ecological environment of this area.