Trends in Forest Greening and Its Spatial Correlation with Bioclimatic and Environmental Factors in the Greater Mekong Subregion from 2001 to 2020

Abstract

Understanding trends of vegetation evolution and its spatial characteristics is critical for sustainable social development in the Greater Mekong Subregion (GMS), which is densely populated and still has uneven economic development. Through Theil–Sen/Mann–Kendall tests, polynomial regression and bivariate local autocorrelation analyses, we investigated vegetation greening trends and their spatial correlation with bioclimatic and environmental variables. The study yielded the following results: (1) Land cover in the GMS has changed significantly over the last 20 years. Conversion between forest and grassland was the main type of change. (2) The upward trend in the forest enhanced vegetation index (EVI) significantly exceeded the downward trend in countries over 20 years. In GMS, the spatial variation in forest trend slope values ranged from −0.0297 a⁻¹ to 0.0152 a⁻¹. (3) Anthropogenic activities have played an important role in forest greening; planted, plantation and oil palm forests exhibit the largest contributions to greening. (4) Changes in forest EVI were most spatially correlated with radiation (12.19% for surface net solar radiation and 12.14% for surface solar radiation downwards) and least spatially correlated with seasonality precipitation (8.33%) and mean annual temperature (8.19%). The results of the analysis of EVI trends in vegetation and their spatial correlation with bioclimatic and environmental variables can provide a reference for strategies aimed for protecting the vegetation ecology.

1. Introduction

Climate change is expected to affect the growth state of vegetation, which will in turn cause regional variations in geographical and environmental elements [1]. This will inevitably affect human productivity and livelihoods, as well as the sustainable development of society. Owing to the spatial heterogeneity of hydrothermal conditions and topography, vegetation responses to climate change show considerable regional differences [2]. Therefore, a long-term analysis of vegetation dynamics and their interrelationships with climatic and geoenvironmental factors at the regional scale is of practical significance importance in ecological construction, carbon balance regulation, and the sustainable use of natural resources in the region.

To date, several vegetation indices, including the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI), have been developed, and various data products have emerged. This includes advanced very-high-resolution radiometer (AVHRR) NDVI, SPOT vegetation (VGT) NDVI, moderate resolution imaging spectroradiometer (MODIS) NDVI/EVI [3], and Landsat NDVI/EVI. They are characterized by a combination of red light absorption and high near-infrared reflection properties for green vegetation leaves and indicate vegetation cover, growth status, biomass, and photosynthetic intensity. They are essential for elucidating the impacts of global change on terrestrial ecosystems and are widely used in regional and global vegetation studies [4,5,6,7]. In contrast with NDVI, EVI is insensitive to the soil background and more sensitive to changes in the dense canopy than NDVI [8,9]. Compared with AVHRR, Landsat, and SPOT-VGT, MODIS EVI has a higher temporal resolution. This allows it to easily overcome the effects of cloud cover and makes it more suitable for analyzing long-term series of vegetation change and responses to climate at the regional scale.

Studying vegetation change in a long time series includes temporal trends in vegetation, the persistence of future changes, mutation analysis, and spatial heterogeneity. Multi-year change trend analysis methods include linear regression [10], Mann–Kendall trend test [11], Sen Slope calculation, and Sen+Mann–Kendall [12]. These have been widely applied in the Yellow River basin and the Beijing–Tianjin–Hebei region of China, as well as in the African and North American continents [13]. However, this approach is limited by its implicit assumption that vegetation changes linearly over a long period of time. The results reflect monotonic changes in vegetation trends, such as monotonic increase and decrease. In this context, it may fail to detect areas with different directions and rates of change [14], or may obscure the presence of short-term “greening” or “browning” patterns through averaging [15]. To explore mutations in vegetation, the segmented linear regression model [16], BFAST [17], and DBEST [18] are popular methods. However, short-term oscillations and abrupt change points may also be caused by orbital drift and instrument calibration, which affect the reliability of the long-term series trend analysis [19]. Furthermore, the impact of climate change on vegetation is mostly a gradual process [19] and can vary nonlinearly with different magnitudes over a range of spatial and temporal scales. This necessitates a detailed analysis of potential changes [14].

The factors affecting vegetation change from climatic constraints, anthropogenic activities, and topography have been investigated [20,21,22,23,24]. Among all factors, temperature and precipitation are the main shaping vegetation change [25,26]. In contrast, bioclimatic variables, which have more ecological importance than mean annual temperature and precipitation, are rarely mentioned [27]. Methods based on multiple regression and residual analysis are widely used to separate and assess the contributions of climatic and anthropogenic factors to vegetation change [28,29,30]. When identifying the magnitude of impacts caused by anthropogenic activities, it is difficult to distinguish between natural and anthropogenic vegetation changes due to climate variability [31,32]. This means that the final contribution of anthropogenic factors obtained is unclear. In addition, vegetation variation and its response to climate change show considerable spatial variations in different regions, owing to different hydrothermal conditions and vegetation types [33]. For example, the vegetation greenness trend in Southeast Asia during the dry season is opposite to that of the Amazon Forest. This highlights the complex effects of climate variability on vegetation dynamics [34]. Therefore, to overcome the shortcomings of linear relationships in regression and residual analyses, a spatial approach is needed to quantitatively assess the correlations between different factors on vegetation change at the regional scale.

The Greater Mekong Subregion (GMS) is located at the junction of Southeast Asia, South Asia, and Southwest China. It is an economic cooperation mechanism initiated by six countries in the Lancang–Mekong River Basin under the initiative of the Asian Development Bank to promote economic and social development in the region [35]. The land cover of the GMS is dominated by forests and agriculture. The rapid development of this region poses a serious threat to its ecosystems. World Wide Fund for Nature (WWF) has stated that this region is vulnerable to the effects of global climate change [36]. This study aims to analyze the trends and spatial patterns of vegetation and land cover in the GMS in the context of regional climate change. The effects of climate and geoenvironmental factors on the spatial and temporal variations of forests in the GMS are also examined. In this study, MOD13Q1 EVI time series data from 2001 to 2020 were obtained from the Google Earth engine. The spatial and temporal variations of vegetation EVI in GMS from 2001 to 2020 were analyzed by combining the monotonic trend and polynomial trend extraction methods. The three main contributions of this study are: (1) extraction and analysis of land cover transformations and EVI trends in the GMS over the past 20 years; (2) comparative analysis of the differences between monotonic and polynomial trend methods; and (3) analysis of the spatial correlations between forest EVI and climate elements as well as geoenvironmental factors.

2. Materials and Methods

2.1. Study Area

The GMS is located in the southwest of China, in Southeast Asia (Figure 1). It was first proposed by the Asian Development Bank in 1992 at a meeting of six countries (regions) in the GMS in Manila, and includes Yunnan Province (CN) of China, Myanmar (MM), Laos (LA), Thailand (TH), Cambodia (KH) and Vietnam (VN). In 2005, the Guangxi province in China joined the GMS as an economic cooperation region. The GMS covers an area of 2,568,600 km², and has a population of approximately 320 million. The terrain is high in the north and low in the south, with numerous mountains and plateaus. The plains in this region are predominantly located in the coastal areas in the southeast, mainly in the alluvial plains and deltas with a wide area downstream of the rivers. The Ayeyarwady, Nujiang–Salween, Lancang–Mekong and Jinsha Rivers are distributed from north to south. The GMS is located in the tropical monsoon region of Asia. From May to October, it is influenced by the southwest monsoon from the sea, resulting in wet and rainy weather. From November to mid-April, it is influenced by the northeast monsoon from the mainland, producing dry and less rainy conditions.

2.2. Dataset

The data used for the study include the EVI, land cover, precipitation, temperature, radiation, topography, and forest management type.

(1) Given the dense vegetation cover in the study area, EVI data were used for vegetation change analysis. The MODIS MOD13Q1 data products had a spatial resolution of 250 m and a temporal resolution of 16 days. They covered the period January 2001 to December 2020, and eight regions of the GMS: h26v06, h26v07, h27v06, h27v07, h27v08, h28v06, h28v07, and h28v08. They were processed and downloaded from Google Earth Engine with the data ID MODIS/061/MOD13Q1.

(2) Land cover data were obtained from MODIS MCD12C1, with a spatial resolution of approximately 0.05°. They were used to analyze land cover change and the EVI trend in the GMS. This land cover dataset comprises 17 major types, including 11 natural vegetation types, 3 land development and mosaic land categories, and 3 non-grassland land type definition categories. They were downloaded from the Google Earth engine.

(3) To analyze the influence of meteorological factors on forest EVI trends, we downloaded meteorological data for temperature, precipitation, and radiation. These data were obtained from the ERA5-LAND reanalysis dataset provided by the European Center for Medium-Range Weather Forecasts (ECWMF). They were obtained using the land-based atmospheric variables simulated by ERA5, and the fifth-generation reanalysis product of ECWMF as forcing factors and simulated using the modified land surface hydrological model HTESSEL and CY45R1 [37]. Compared with ERA5, ERA5-LAND has a spatial resolution of 0.1° (9 km) and a temporal resolution of 1 h [37].

(4) To analyze the spatial correlations between environmental factors and EVI trends, we downloaded ASTER GDEM V2 [38], which was calculated based on data from the Advanced Spaceborne Thermal Emission and Anti-Emission Instrument (ASTER). ASTER GDEM V1 was released on 29 June 2009. The GDEM V2 solved the anomalies in the GDEM V1 data and improved the accuracy of spatial resolution and elevation.

(5) To differentiate anthropogenic activities in the analysis of the forest EVI trends, we downloaded globally consistent forest management map data (FMM) [39]. This map was derived from the 226 K reference dataset from Geo-Wiki with 2015 PROBA-V satellite imagery at a resolution of 100 m. It has an accuracy ranging between 58% and 80% for the forest management classes. This map has a high level of spatial detail for the most prevalent forest management classes, such as naturally regenerating forests without any signs of management (NRF), naturally regenerating forests with signs of forest management (NRFM), planted forests (PF₁), plantation forests (PF₂) with a rotation of up to 15 years, oil palm plantations (OPP), and agroforestry (AF).

2.3. Methods

The flow of the research method is illustrated in Figure 2. After the performing pre-processing steps, such as EVI filtering, maximum value extraction, land cover reclassification, and calculation of biological variables and annual mean of net surface solar radiation, the monotonic and polynomial trends were extracted from the EVI series data. The spatial statistics were then determined. Finally, bivariate local Moran’s I (BiLISA) was used to analyze the spatial correlations between the climate and the environmental variables and EVI trends. The main steps are as follows.

(1)

Remote sensing data pre-processing

First, Savitzky Golay (S–G) filtering was performed on MODIS EVI series data from 2001 to 2020 using the Google Earth engine to remove outliers. The de-clouding of MODIS EVI was based on the labels of MODIS clouds. The maximum EVI value of each pixel was extracted from the valid EVI series data using the maximum value synthesis method [40] to obtain the $EV{I}_{max}$ in the ith year. The $EV{I}_{max}$ represented the vegetation cover and served to eliminate the influence of clouds and water vapor on the EVI and improve the reliability of the data.

In this study, the original 17 classes of MCD12Q1 were reclassified into seven classes: water bodies, forests, grasslands, croplands, cropland/natural vegetation, urban, and bare land. Evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, mixed forests, and closed shrublands were classified as forests. Open shrublands, woody savannas, savannas, grasslands, and permanent wetlands were classified as grasslands. Permanent snow/ice and water bodies were classified as water bodies (Appendix A,

Table A1. MODIS land cover reclassification.

Classes	Code	Reclassified Classes	New Code
Evergreen needleleaf forests	1	Forests	1
Evergreen broadleaf forests	2
Deciduous needleleaf forests	3
Deciduous broadleaf forests	4
Mixed forests	5
Closed shrublands	6
Open shrublands	7	Grasslands	2
Woody savannas	8
Savannas	9
Grasslands	10
Permanent wetlands	11
Croplands	12	Croplands	3
Urban and built-up lands	13	Urban and Built-up Lands	4
Cropland/natural vegetation mosaics	14	Cropland/Natural Vegetation Mosaics	5
Barren	15	Barren	6
Permanent snow and ice	16	Water Bodies	7
Water bodies	17
Unclassified	255	Unclassified

).

For the ERA5-land hourly data, nine bioclimatic variables were calculated from the ERA5 land hourly data according to the definition of the “biovars” function in the “dismo” package in R 4.2.1. We then calculated the average values of five climatic variables, namely surface net solar radiation (SSR), surface solar radiation downward (SSRD), soil temperature layer 1 (STL1), volumetric soil water layer 1 (SWVl), and total evaporation (E) for 2001–2020. These data were then re-projected, cropped, and resampled into the spatial projection and resolution of MODIS EVI for subsequent data analysis. Finally, three topographic variables (elevation, slope, and aspect) were calculated using GDEMV2. These data are listed in Appendix A,

Table A2. Biometeorological and topographic factors.

Variable	Code	Unit
Annual mean temperature	bio01	°C/a
Temperature seasonality	bio04	°C/a
Max Temperature of warmest month	bio05	°C/a
Min temperature of the coldest month	bio06	°C/a
Mean temperature of wettest quarter	bio08	°C/a
Mean temperature of driest quarter	bio09	°C/a
Mean temperature of warmest quarter	bio10	°C/a
Annual precipitation	bio12	mm/a
Precipitation seasonality	bio15	%/a
Elevation	elevation	m
Topographic Slope	slope	°C
Surface net solar radiation	ssr	J m−2
Surface solar radiation downwards	ssrd	J m−2
Soil temperature level 1	stl1	K
Volumetric soil water layer 1	swvl1	m3
Total evaporation	e	m of water equivalent

.

(2)

Analysis of the EVI trend of vegetation

In this study, the Theil–Sen + Mann–Kendall (TS–MK) [41] trend analysis method was used. The TS–MK method uses Theil–Sen to calculate the trend slope and Mann–Kendall to test the significance of long time series trends. Mann–Kendall is a nonparametric statistical method widely used for the trend discrimination of meteorological and hydrological data. The Mann–Kendall method has no mandatory distribution requirements for series data and is not limited by anomalies. It is performed to statistically assess whether there is a monotonically increasing (upward) or decreasing (downward) trend in the EVI over time. A monotonically increasing trend implies that the variable increases over time. However, this trend may or may not be linear.

The Theil–Sen trend slope is calculated as follows:

$$\beta =Median\left(\frac{x_j-x_i}{j-i}\right) \forall j>i$$

(1)

where $1<i<j<n$, $\beta$ is the slope, $\beta >0$ indicates an upward trend in vegetation, and $\beta <0$ indicates a downward trend in vegetation.

The statistic $S$ for the Mann–Kendall trend test is calculated as follows:

$$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}sgn\left(x_j-x_i\right)$$

(2)

where

$$sgn\left(x_j-x_i\right)=\{\begin{array}{c}1 x_j-x_i>0\\ 0 x_j-x_i=0\\ -1 x_j-x_i<0 \end{array}$$

(3)

The length of the time series in this study is $n$ = 20 (2001–2020), and the statistic $S$ approximately obeys the standard normal distribution. Therefore, statistic $Z$ is used to perform the trend test, and the significance level is considered to be $\mathsf{\alpha }=0.05$, $Z_{1-\alpha }=Z_{0.975}=1.96$.

The statistic Z is calculated as follows:

$$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}$$

(4)

where

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

(5)

where $n$ is the number of data points in the sequence, $m$ is the number of knots (groups of recurring data points) in the sequence, $t_i$ is the width of the knot, that is, the number of repeated data points in the $i$th group of repeated datasets. A bilateral trend test was performed for the $Z$ values. $\left|Z\right|\ge 1.96$ indicates that the trend of the time series passed the significance level of 0.05.

Considering that TS–MK only detects monotonic trends, we used a polynomial trend [14] to detect nonlinear features in the EVI series data. The polynomial trend is a three-stage procedure for generating linear, quadratic, cubic, concealed, and non-trending classes of vegetation index time series. Polynomial trends are extracted using a recursive judgment method by statistically testing the significance of linear, quadratic, and cubic trends. During the first two significance tests, pixels were assigned cubic or quadratic trend classes when passing the test, but were treated as concealed classes if the slope coefficients were not statistically significant. Finally, the linear trend class was assigned to a pixel in the case of passing the significance test. Otherwise, it was regarded as having no trend.

(3)

Spatial correlation analysis of EVI trends with climate and environmental variables

To detect the spatial correlation characteristics of trends in forest EVI with bioclimatic and environmental variables, we used the BiLISA tool in GeoDa. BiLISA can reveal the spatial clustering and spatial differentiation characteristics of one attribute in each spatial cell and another attribute in the neighboring cells. BiLISA was calculated as follows:

$$I_{pq}=W_p{\displaystyle \sum }_{q=1}^n{Z}_{pq}{W}_q$$

(6)

where $I_{pq}$ is the local spatial correlation index in spatial cells $p$ and $q$, and $W_p$ and $W_q$ are the normalized values of the variance of the attribute values in spatial cells p and q, respectively. BiLISA can generate five types of spatial associations: spatially nonsignificant, high–high clustering (HH), high–low clustering (HL), low–low clustering (LL), and low–high clustering (LH). A total of nine biometeorological, two topographic, and five environmental variables were used in this study. To explore the correlation of these variables with EVI trends, we calculated their slopes for 2001–2020 using linear regression.

3. Results

3.1. Characteristics of Land Use/Land Cover Changes in the Greater Mekong Subregion

The spatial distribution of land cover of the GMS during 2001–2020 and its changes are shown in Figure 3. During 2001–2020, the land cover changes in the GMS were relatively significant, and 81.88% of the land cover types remained spatially stable. Forests and grasslands were the dominant land cover types in the GMS (Figure 3), accounting for more than 78% of the entire area. Among the categories where land cover change occurred, the pronounced ones were the conversion from forest to grassland (171,656.09 km²) and conversion from grassland to forest (114,631.71 km²). Conversion from forest to grassland was mainly distributed in Laos and Myanmar. Meanwhile, conversion from grassland to forest was mainly distributed in Guangxi and Yunnan, China, and the northern part of Thailand. In terms of forest loss and gain, Guangxi and Yunnan in China had the largest forest gain (64,036.84 km²), while Laos (3786.70 km²) and Cambodia (2695.88 km²) had the smallest gain. The largest forest loss was in Myanmar (48,560.93 km²) and the smallest was in China (14,509.18 km²). In terms of grassland loss and gain, Myanmar had the largest grassland gain (54,315.94 km²), and the remaining countries exhibited a small differences in grassland gain (33,511.51–37,695.95 km²). China had the largest grassland loss (78,044.89 km²), of which 49,281.51 km² was transformed into forest, and the smallest grassland losses were in Laos (6636.01 km²) and Cambodia (11,674.21 km²).

3.2. Monotonic Trends in EVI of Vegetation with No Change

From the results of the MK significance trend test, the non-significant pixels encompassed a larger proportion of the area of each country, and were ranked according to value: Laos (89.7%) > Myanmar (83.3%) > Cambodia (74.1%) > Vietnam (73.2%) > Thailand (71.3%) > China (66.4%). Among the significant pixels, the largest proportion of the downward trend was Vietnam (5.9%), and the smallest proportion was Thailand (1.8%). In terms of the proportion of the area in an upward trend, Yunnan and Guangxi provinces in China had the largest increase (31.5%) and Laos the smallest (9.5%). Spatially (Appendix A, Figure A1), the areas with significant downward trends were mainly located around cities, such as the Red River Delta in Vietnam, Phnom Penh in Cambodia, Bangkok in Thailand, and Nanning City in Guangxi, China. The upward trend was concentrated and continuous in Guangxi Province, China, the Lower Mekong Plain, the Lower Ayeyarwady River in Myanmar, and northeastern Thailand.

Figure 4a shows the distribution of significant trends in forest EVI without type change. Areas with a significant increase in forest EVI were located in several mountains, including the Truong Son Mountains, Birao Mountains, Arakan Mountains, and Shan Plateau as well as eastern Guangxi and southwestern Yunnan, China. Areas with a significant downward trend accounted for 8.1% of the area of the upward region and were spatially dispersed. There were also areas of intense downward movement in certain areas, such as Khammoung and Xieng Khouang in Laos. Areas with significant increases in EVI for grasslands without type change were mainly located in Guangxi and western and northern Yunnan in China, and northeastern and southern Thailand. A downward trend is present for in Laos, western Vietnam, and northern Cambodia.

Figure 5 shows the distribution of the area and slope of the significant trend of EVI for unchanged forests and grasslands. Figure 5a shows that the largest area of significant EVI increase was in China, followed by Myanmar, and the smallest was in Cambodia. The upward trend for the EVI in the forests in each country was greater than the downward trend. For grasslands (Figure 5b), the largest area of significant EVI increase was in China and the lowest in Laos. Figure 5c,d shows the slope distribution of areas with significant EVI trends for forests and grasslands. In terms of forests, spatial variation in the slope values of forests in the GMS ranged from −0.0297 a⁻¹ to 0.0152 a⁻¹ over the last 20 years. For the spatial distribution, the largest upward slope was observed in China with a mean value of 0.0077 a⁻¹, and the smallest was observed in Vietnam with a mean value of 0.0058 a⁻¹. Dispersion of the upward slope was strongest in China and Cambodia. Meanwhile, the mean of the downward slope was more variable in the forests and more discrete than the upward trend. In terms of grasslands, the largest upward slope was for Cambodia (0.0081 a⁻¹) and the smallest was for Laos (0.0062 a⁻¹). The dispersion of the downward slope was greater than that of the upward slope. In summary, the overall trend for forests and vegetation without type change in the last 20 years in the GMS was upward.

3.3. Polynomial Trends in EVI of Vegetation with No Change

We extracted polynomial trends for forests and grasslands that had not experienced type changes (Figure 6). Five types of results were included in the polynomial trend algorithm, namely linear, quadratic, cubic, concealed trend, and no-trend classes. Only the first four types are presented here. Figure 6a shows the polynomial trends of the forests over the last 20 years. Spatially, the polynomial trends for forests were concentrated in the Shan Plateau, Arakan Mountains, Birao Mountains, Xieng Khouang Plateau, and lower and middle Mekong Mountains. Concealed forest trends were widely distributed and dispersed. Linear trends in forests were mainly observed in Ha Tinh Province and the Birao Mountains in Vietnam, the Cardamom Mountains in Cambodia, and the eastern and southern parts of Guangxi, China. Meanwhile, dispersed linear trends were observed in the Shan Plateau in northern Myanmar and the Hengduan Mountains in Yunnan. Quadratic and cubic trends were not spatially significant. Figure 6c shows that the forest had the largest proportion of concealed trends in terms of area, followed by linear trends, while the proportion of quadratic and cubic trends was smaller. Among the six countries, the concealed and linear trends were the largest in Myanmar. In China, the largest proportion of linear trends were found in forests, followed by concealed trends. Cambodia had the smallest area for all trends. In terms of grasslands, linear trends were mainly distributed in the eastern part of the Korat Plateau and the southern part of Changwat Krabi in Thailand, Guangxi province and the eastern part of Yunnan province in China, and Mon State in Myanmar (Figure 6b). The concealed, quadratic, and cubic trends were spatially dispersed, with only the quadratic trend being more concentrated in the northern part of Ho Chi Minh City in Vietnam. Figure 6d shows the area proportion for each trend. Compared with forests, the linear trend was larger in grasslands, followed by concealed trends. There was little difference in the proportions of quadratic and cubic trends for grasslands among the six countries.

Compared to TS–MK (Figure 5), the polynomial trend method (Figure 6) can detect more pixels with trends. In terms of forests, 11.28% of the pixels that did not show a significant trend in TS–MK were judged as concealed trend classes by polynomial trends, and 3.02% of the pixels were judged as linear (2.49%), quadratic (0.25%), and cubic (0.28%) trend classes. In pixels with TS–MK upward trends, the percentages of concealed, linear, quadratic, and cubic trend classes were 2.07%, 68.56%, 4.48%, and 4.8%, respectively. For the pixels with TS–MK downward trends, the percentages of concealed, linear, quadratic, and cubic trend classes were 2.92%, 58.69%, 6.4%, and 8.7%, respectively. With grasslands, the polynomial trend determined the pixels in which the TS–MK method showed no significant trend as concealed (14.74%), linear (3.5%), quadratic (0.62%), and cubic (0.49%) trend classes. The percentage of linear, quadratic, and cubic trend classes was higher than that of the forest in pixels with TS–MK upward and downward trends, whereas the percentage of concealed trends was lower than that of forests.

3.4. Impact of Human Activities on EVI Trends in Forests

Figure 7a shows the distribution of significant trends in EVI for the forest management types defined by the FMM and the proportion of their area in each category. The areas of NRF and NRFM areas were the largest in the GMS. The proportion of significant trends was the largest in PF₂ and PF₁. This indicates that anthropogenic activities in the GMS strongly influence significant changes in forest EVI. For each country, Figure 7b,c shows the proportions of upward and downward trends in EVI for each forest management class. In terms of upward trends (Figure 7b), PF₁ and PF₂ had the largest proportions in China, reaching 27.76% and 40.46%, respectively, whereas AF had the lowest (9.53%). In Vietnam, the OPP had the largest upward trend, whereas the AP had the smallest. In Thailand, PF₂ showed the largest upward trend and NFRM the smallest. In Myanmar, the OPP showed the largest upward trend and the AP the smallest. In Cambodia, PF₂ showed the largest upward trend and NFRM the smallest. Overall, China showed relatively high upward trends in NRF, NRFM, PF₁, and PF₂, whereas other countries showed high upward trends in PF₂ and OPP. As shown in Figure 7c, Vietnam showed the most pronounced downward trends in PF₂ (2.11%), OPP (3.46%), and AF (4.22%), while other countries showed the largest proportion of downward trends in AF. However, the downward trend in GMS forests affected a small proportion (ranging from 0.36% to 4.2%). Figure 7d,e show the slope of the significant trend in EVI for different forest management types. In the upward trends, the PF₂ of GMS showed the largest variation ($0.0075\pm 0.0027$), and the smallest was AF ($0.0060\pm 0.0023$). For the downward trends, PF₁ had the highest variation ($-0.0071\pm 0.0044$), and its dispersion was larger. PF₂ exhibited the smallest variation ($-0.0055\pm 0.0023$).

3.5. Spatial Correlation Analysis of Trends in EVI with Climate and Environmental Variables

Overall, among the 16 variables, the average percentage of those without significant clustering characteristics with forest EVI trends reached 60.37%. Meanwhile, among those with significant clustering characteristics with forest EVI trends (p < 0.05), the average percentages of high–high clustering (HH), low–low clustering (LL), low–high clustering (LH), and high–low clustering (HL) were 9.56%, 9.86%, 11.02%, and 9.19%, respectively. In terms of clustering characteristics, the largest percentage of EVI with no significant characteristics was observed for bio01 (67.25%) and the smallest for SSR (51.24%). In HH clustering, the largest proportion was observed for EVI with a topographic slope (15.15%) and the smallest in bio12 (3.47%). In LL clustering, the largest proportion was observed for EVI with elevation (15.61%) and the smallest for topographic slope (1.50%). In LH clustering, the largest proportion was observed for EVI with a topographic slope (16.29%) and the smallest with bio12 (6.93%). In the HL clustering, the largest proportion was obtained for EVI with elevation (13.59%) and the smallest with topographic slope (1.19%). Spatially, LL clustering was most evident in the northern part of Myanmar, as shown in Figure 8e,g–j,o. LL clustering also exists in transnational regions, such as the borders of China, Vietnam and Laos, as shown in Figure 8c,h,i,p, and the borders of Vietnam, Laos and Cambodia, as shown in Figure 8b,j,l,m. LH clustering is mainly present in the northern part of Myanmar, as shown in Figure 8a–g,p. HL clustering was most evident in Yunnan and Guangxi, China, and Myanmar, as shown in Figure 8a,c–f,h–n. HH clustering was mainly distributed in northern Myanmar and southern Yunnan, China (Figure 8a,d,e,g,j–m,p), central Truong Son Mountains of Laos and Vietnam (Figure 8g–i,k), and the Birao Mountains of Myanmar and southern Thailand (Figure 8b,d–f,n,o).

4. Discussion

Overall, our results showed (Figure 3) that forests were decreasing in five of the six GMS countries, with Laos and Cambodia showing the largest reductions over the last 20 years. Meanwhile, China was the only country with an increase in the forest area. According to Tölle [42], forest dynamics occurred in stages, with a sharp decrease in forest area after 2000, which continued until 2018, as in Cambodia and Vietnam. In recent years, some countries have shown signs of reforestation or slow deforestation, such as Thailand [43,44]. According to food and agriculture organization studies, the recovery of forest cover in some countries is attributed to the expansion of plantations [45,46]. Therefore, despite the general slowdown in deforestation and efforts to increase forests in these countries, the trend of tropical forest decline has not changed [47]. Much of this forest decline has shifted to grasslands, as seen in southern Cambodia and northern Laos (Figure 3). These declines occurred in the marginal parts of forested areas, which is consistent with the trend observed by Tölle [42]. Land policies in China have been adjusted in recent years. In Guangxi, during the Tenth Five-Year Plan, forestry projects included the return of cultivated land to forest, the construction of fast-growing and productive forest bases, and the construction of nature reserves, Pearl River protection forests, and coastal protection forests [48]. Open forest land with 10–30% depression and high-cover grassland with >50% cover were planted into timber and economic forests [49]. This is the main anthropogenic factor affecting forest growth in China.

In recent years, global vegetation is becoming greener [50], including in India, southern and southeastern China, and Southeast Asia [51]. Study timespans are important for determining time-series trends in vegetation indices [16]. It has been shown that global vegetation trends after 2000 are highly uncertain [47]. Multiple trend detection methods are required to assess the coherence of trends. In this study, a monotonic trend assessment method and polynomial trend analysis were used. With TS–MK, the polynomial trend was able to find more trend pixels in EVI, especially those with concealed trends. In contrast, 58.69–68.56% of the significant trend pixels were linear. This does not mean that TS–MK overestimated the trend variation of the EVI. The TS–MK method is insensitive to noise, does not require the data to obey a normal distribution, and has high reliability in monotonic trend detection. The difference in the number of linear trends may be attributed to the polynomial trend that assigns pixels from the monotonic trend to the quadratic and cubic trend classes. In terms of spatial distribution, the trends were similar. In the case of forests, areas where trends were clearly concentrated were distributed in eastern China and central Laos. Meanwhile, areas where trends appeared more dispersed were distributed in the mountainous areas between Myanmar and Yunnan, China. The similarity in spatial distribution characteristics of the results of the two detection methods was more evident in grassland trend maps. From the results of vegetation index trends, all six countries located in the study area had unique land use policies, and therefore, may have unique vegetation trend trajectories. Chen et al. [47] showed that the “greening” of Earth is mainly attributed to increasing in human agriculture and forest replanting. In this study, we excluded agriculture and areas where forest categories changed over a 20-year period. As shown in Figure 5 and Figure 6, these increasing trends in EVI were significantly greater than the decreasing trends. However, spatially, the increasing and decreasing trends in EVI are intertwined, which is consistent with the observations of Chen [51].

In terms of land cover change, the GMS had the greatest reduction in forests, and the greatest increase in grasslands (Figure 3), and a large proportion of conversion from forests to grasslands. One of the main driving factors is deforestation by agricultural production in Southeast Asian countries [52]. The high value of forests to the ecology of the GMS and the sustainability of human society is well known. Therefore, this study used BiLISA to analyze the spatial clustering characteristics of biometeorological and topographic factors and forest EVI trends. From the results, 16 variables with no significant spatial correlation with forest EVI trends in the GMS accounted for most of this region. Net surface solar radiation (SSR; 12.19%), downward surface solar radiation, (SSRD; 12.14%), and elevation (11.29%) were the three variables with the largest average distribution correlated with EVI trends (HH, LL, LH, and HL). This was followed by STL1 (soil temperature class 1, 10.52%), Bio05 (maximum temperature in the warmest month, 10.47%), Bio09 (average temperature in the driest quarter, 10.36%), and Bio08 (average temperature in the wettest quarter, 10.26%) in the biometeorological variables, although the differences in the average proportions between them were relatively small. The smallest average proportions of spatial clustering characteristics with EVI trends were the topographic slope (8.53%), Bio15 (precipitation seasonality, 8.33%), and Bio01 (annual mean temperature, 8.19%). Nemani [20] suggested that vegetation growth in Southeast Asia is mainly limited by radiation rather than by other climatic factors such as precipitation and temperature. Results of our BiLISA analysis also showed that radiation (SSR and SSRD) had the largest spatial correlation with forest EVI trends. In addition, forest EVI trends were also influenced by elevation due to differences in hydrothermal, light, and atmospheric conditions resulting from different elevations. Our results further confirmed that the precipitation and temperature variables Bio15 and Bio01 were the two variables with the lowest spatial correlation with forest EVI trends; however, seasonal changes in precipitation and temperature were more spatially correlated with forest EVI trends.

To date, we have only analyzed monotonic and polynomial trends, which complicated determining of the response of vegetation to weather extremes. In future studies, we will use long-term time-series datasets to capture trends and their turning points as well as related influencing factors. In this study, we analyzed anthropogenic activities in unchanged forest areas. Considering the higher anthropogenic activities in the GMS, the impacts of anthropogenic activities on other vegetation types will be further explored.

5. Conclusions

This study aimed to analyze the change characteristics of land cover as well as grassland and forest EVI of six countries in the GMS, and their spatial correlation with bioclimatic and environmental variables. We used TS–MK and polynomial regression to extract the vegetation trends in MODIS EVI over the past 20 years, and BiLISA was used to characterize the spatial clustering of bioclimatic and environmental variables with forest EVI trends. The results of the study show that land cover changes in the GMS were relatively large during 2001–2020. Among types of land cover change, conversion between forest and grassland was the main trend. The EVI trend in forests and grasslands without a type change in the GMS was upward. Forest EVI trends in the six countries were strongly influenced by anthropogenic activities, especially planted forests and plantation forests in China, and oil palm plantations in the other countries. Changes in forest EVI were most spatially correlated with radiation and least spatially correlated with the seasonality of precipitation and mean annual temperature. Considering the relatively high population density and unbalanced economic development in the GMS, the sustainable development of society is urgently needed for soil and water conservation, carbon sequestration capacity, and protecting the economic value of forests. While working towards, regional economic development, countries should work together to strengthen the regulation of anthropological activities and effectively reduce forest loss.