Characterizing the Long-Term Landscape Dynamics of a Typical Cloudy Mountainous Area in Northwest Yunnan, China

Abstract

Detailed knowledge of landscape dynamics is crucial for many applications, from resource management to ecosystem service assessments. However, identifying the spatial distribution of the landscape using optical remote sensing techniques is difficult in mountainous areas, primarily due to cloud cover and topographic relief. Our study uses stable classification samples from mountainous areas to investigate an integrated approach that addresses large volumes of cloud-cover data (with associated data gaps) and extracts landscape time series (LTS) with a high time–frequency resolution. We applied this approach to map LTS in a typical cloudy mountainous area (Erhai watershed in northwestern Yunnan, China) using dense Landsat stacks, and then we also used the classified results to investigate the spatial–temporal landscape changes in the study area at biennial intervals. The overall accuracy of the landscape classification ranged from 81.75% to 88.18%. The results showed highly dynamic processes in the landscape throughout the study period. Forest was the main land cover type, covering approximately 39.19% to 41.68% of the total study area. Alpine meadow showed fluctuating trends, with a net loss of 11.22% and an annual reduction rate of −0.4%. Shrub cover increased by 1.26%, and water bodies showed a small decrease in area, resulting in an overall net change of −0.03%. Built-up land and farmland areas continued to expand, and their annual growth rates were 1.52% and 1.06%, respectively. Bare land showed the highest loss, with a net change of 228.97 km². In the Erhai watershed, all the landscape classes changed or transitioned into other classes, and a substantial decrease in bare land occurred. The biennial LTS maps allow us to fully understand the spatially and temporally complex change processes occurring in landscape classes; these changes would not be observable at coarse temporal intervals (e.g., 5–10 years). Our study highlights the importance of increasing the temporal resolution in landscape change studies to support sustainable land resource management strategies and integrate landscape planning for environmental conservation.

1. Introduction

A landscape is a complex system comprising a geographic and/or an ecological space in which organisms live. Landscapes are affected by external disturbances (such as earthquakes, erosion, climate change, fire, human disturbances, etc.) and internal perturbations (such as physical flow of matter and energy, community turnover, etc.) with change processes that evolve over time. Landscape dynamics are very important for environmental management and conservation; they are a vital field of sustainable development research [1]. Mountain ecosystems around the world are particularly vulnerable, especially in the plateau mountainous regions of Southwest China, due to climate change and human disturbances [2]. Therefore, a detailed and precise characterization of landscape dynamics is essential for understanding the impacts on the mountainous environment and ecological systems.

As a typical cloudy mountainous region in China, the Erhai watershed (EHW) in northwestern Yunnan is located in a low-latitude, high-altitude belt [3,4]. In this area, complex topography forms steep environmental gradients that generate climatic differences and landscape diversity. The EHW is a typical hybrid zone, with a composition of lakes, mountains, cities, and numerous types of vegetation, and the landscapes in this region have a vertical zonal distribution. The EHW is recognized as a global biodiversity hot spot and supports rich wild fauna and flora [4,5,6]. Due to its unique geographical environment, the EHW is especially vulnerable to climate change and anthropogenic disturbances; climate factors exert strong control on mountain landscapes, and anthropogenic behavior impacts the landscape dynamics of intermontane watershed areas. Hence, terrestrial and aquatic ecosystems in the EHW are threatened by climate change and anthropogenic-induced changes in landscape dynamics [5,7], which lead to the alteration of ecological processes and ecosystem services [4], habitat degradation or loss, and landscape fragmentation [7]. Knowing when and where changes in the landscape of the EHW have occurred is fundamental to assessing land management impacts and to better understanding how the region’s ecosystem responds to climatic variability and anthropogenic pressures; this knowledge can answer scientific and ecological questions. Therefore, the high temporal resolution landscape time-series (LTS) maps over a long-term period are urgently needed to analyze how the landscape evolves in the EHW.

Satellite remote sensing can monitor landscape dynamics and environmental conditions at broad spatial scales and high temporal frequencies over long periods [8,9,10]. However, the use of optical satellite-based imagery (such as Landsat imagery) is onerous in mountainous areas due to persistent heavy cloud cover, data gaps, and topographic variations [10,11]. Thus, the development of LTS maps to analyze landscape dynamics from optical satellite-based images faces numerous challenges: (1) consecutive cloud cover may obscure information related to seasonal or inter-yearly landscape change [12]; (2) topographic influences are obvious in mountainous areas, and impose additional differences in spectral radiance within any particular land cover due to changes in surface slope angle and aspect, leading to errors that may impact landscape classification and pattern analysis [13]; and (3) little or no historical ground truth or other reference data for landscape classification samples (training or validation samples) are available [14]. Consequently, determining LTS changes across a long time span is a difficult task in cloudy mountainous areas. Various approaches have been proposed to reconstruct satellite data to minimize the impact of cloud cover and topographic relief, such as pixel-based compositing [15,16], spatial–temporal gap filling [17,18], and temporally composited mosaicking [19]. These studies provide many reference solutions to address data gaps. However, the lack of historical ground truth or other reference data is still the main obstacle to mapping long-term LTS with a high temporal resolution (annual or biennial interval).

In view of the urgent need for information on landscape change, this study mapped high temporal resolution LTS as a starting point to analyze and understand the long-term landscape dynamics that occurred in the EHW. For this purpose, two specific aims were deployed to analyze the long-term spatiotemporal dynamics of the landscape in the study area:

• (1). To propose a composite approach to address large volumes of cloud-cover satellite imagery and map high time–frequency LTS in mountainous areas.

• (2). To characterize and estimate landscape dynamics using the LTS maps in the EHW.

2. Materials and Methods

2.1. Study Area

The EHW is located in northwestern Yunnan, China, between 25°25′33.891″–26°25′50.961″ N and 99°50′5.854″–100°26′11.627″ E (Figure 1). The EHW area is 2906.36 km², with a maximum altitude of 4050 m, a minimum altitude of 1524 m, and a mean altitude of 2776.36 m. The EHW has a subtropical monsoon climate, with a dry season from December to April and a rainy season from May to November.

The study area includes the city of Dali and parts of the city Eryuan in Dali Bai Autonomous Prefecture, Yunnan Province. The EHW is situated in the watershed belt of the Jinsha River, Red River, and Lancang River [3] and is surrounded by mountain chains (the Hengduan range). The natural landscape is the dominant background throughout the EHW, and manual or semiartificial landscapes are embedded in the area [5]. More than 2000 km² of the terrestrial surface is mountainous, and the lake area accounts for more than 8.82% of its total area. Erhai Lake (one of the lakes inside this area), located in the central part of the EHW, is the second-largest freshwater lake [20], and is known as a “pearl” of the Yunnan–Guizhou Plateau. The western side of Erhai Lake is adjacent to the steep Cang Mountains, which rise to more than 2000 m above the lake surface [21]. The EHW is subjected to a subtropical plateau monsoon climate [3], and the climate and landscapes are vertically distributed. The EHW is considered a special area with tremendous biodiversity and landscape diversity.

Dali and parts of Eryuan are situated inside the EHW, and the area surrounding Erhai Lake is one of the most developed regions in Yunnan Province, with a high population density [22]. Moreover, the EHW is well known for its ethnic customs and natural scenery and has become one of most popular tourist destinations in China; it also contains part of the Cangshan Global Geopark.

2.2. Data Used

Since the early 1970s, consecutive Landsat missions have provided remote sensing data with fine spatial and spectral details that have captured long-term landscape characteristics with high time–frequency [23]. However, the revisit period of Landsat satellites (16 days) and the meteorological conditions of the study region (significant cloud cover on many days) reduce the availability of Landsat images. We acquired the cloud-cover statistics automatically calculated by the United States Geological Survey. Figure 2 shows the mean cloud cover per month for the study area from 1986 to 2017. Between June and September, the cloud cover is greatest, and cloud-free or less cloud cover occurs on very few occasions. These states bring many challenges for mapping vegetation landscapes in the leaf-on season because they increase the chance of obtaining numerous cloud-covered Landsat images. Therefore, the images were collected between 1 March and 30 October to represent the growing season. To minimize the influence of clouds and cloud shadows, the maximum cloud coverage was set to 80% [24]. A total of 395 images acquired from 1986–2017 were included in this study (Figure 3). All the imagery data were obtained from the EROS Science Processing Architecture (ESPA) (https://espa.cr.usgs.gov (accessed on 1 July 2019)), including surface reflectance (SR) data, top of atmosphere (TOA) reflectance data, and brightness temperature (BT) data.

To facilitate mapping land cover and to characterize landscape dynamics, a set of auxiliary geospatial datasets were collected for this study, which included: (1) Google Earth maps (2009 to 2017) for validating the accuracy of the stable samples and (2) digital elevation model (DEM) from the Advanced Spaceborne Thermal Emission and Reflection Radiometer.

2.3. Methods

In a mountainous area with cloudy and complex terrain, it is quite challenge to map long-term LTS with a high temporal resolution (annual or biennial interval). This is mainly because there is little ground truth or other reference data for classification samples (training and validation samples). Therefore, a specific technological framework is required to adequately map and monitor landscape dynamics in the EHW over a long-term period (Figure 4). We propose a five-stage framework for mapping biennial interval LTS in the EHW over the period from 1986 to 2017: (1) topographic radiometric correction was conducted to minimize the effect of terrain variation; (2) clouds and shadows were masked out to leave unobstructed areas; (3) feature indices (such as the normalized difference built-up land index (NDBI), normalized difference vegetation index (NDVI), etc.) and maximum value composite (MVC) methods were extracted to fill the missing data pixels; (4) water bodies were masked out, and a set of stable samples were selected (for training and validation) to eliminate the trouble of lacking historical classification samples; and (5) landscape classification and dynamic analysis were conducted.

2.3.1. Topographic Radiometric Correction

All required images from the ESPA were L1T and projected in the UTM coordinate system (WGS 84 datum, Zone 47 North). However, in the study area with strong topographic variations, the classification results produced from the ESPA Landsat data may not be reliable enough. Topographic radiometric correction is a necessary preprocessing step to minimize the effect of topography on surface reflectance. Therefore, a preprocessing step of topographic radiometric corrections was performed on the Landsat data using the SCS+C model to eliminate or reduce the topographic impact on the SR data [13,25].

2.3.2. Cloud Mask

We downloaded images with less than 80% cloud cover for 1 March to 30 October during the 1986–2017 period. Cloud masking is a key process to eliminate clouds and cloud shadows; we employed an algorithm for mountainous areas called Mfmask (mountainous Fmask) that can automatically detect clouds and cloud shadows. This algorithm is built upon the Fmask algorithm and rectifies some limitations of the Fmask algorithm in mountainous areas [26]. TOA reflectance, BT data, and DEM data were used as inputs to the Mfmask algorithm. TOA and BT were directly obtained from the ESPA.

2.3.3. Feature Index Derivation and Image Composition

In this study, landscapes are classified into seven types, namely, forest, alpine meadow, shrub, farmland, water body, built-up land, and bare land. A large number of remote sensing indices for landscape classification and feature extraction have been described [14]. We calculated 12 remote sensing indices to quantify the land surface information. These multispectral indices were chosen to reflect a difference in surface status, such as vegetation features, bare land, built-up land, water bodies, etc. Among them, the vegetation indices include the ratio vegetation index (RVI) [27], NDVI [28], woodiness index (WI) [29], and spectral variability vegetation index (SVVI) [30]; the water indices include the modified normalized difference water index (MNDWI) [31] and normalized difference moisture index (NDMI) [32]; the built-up land indices include the enhanced built-up and bareness index (EBBI) [33], NDBI [34], and urban index (UI) [14]; and the bare land indices include the normalized difference bareness index (NDBAI) [35], normalized difference bare land index (NBLI) [36], and normalized difference soil index (NDSI) [37].

Then, we used the Mfmask layer to mask out all the nonvalid data (clouds and cloud shadows) to calculate various feature index imagery. To minimize the influence of clouds and cloud shadows, the feature index imagery that masked out nonvalid data between 1986 and 2017 was combined using the MVC method at two-year intervals. The biennial filled images of various feature indices were finally obtained in the growth season (1 March and 30 October).

To reduce redundant index information and improve the discrimination of indices for landscape characteristics, samples of typical landscapes were selected from 12 index images, and the capacity of these 12 indices was comparatively analyzed using these samples (Figure 5). In the vegetation index category (RVI, NDVI, WI, and SVVI), vegetation features are enhanced in both the NDVI and RVI images, and alpine meadow and bare land have similar values in these two indices. The WI and SVVI indices can distinguish alpine meadows from bare land, and both indices are differentiated for different landscape types. Two water indices can enhance water information, while others are suppressed. Compared to the NDMI images, the water bodies in the MNDWI images have higher values than the other types. In the NDMI, various landscape types can be distinguished from each other due to different surface humidities. In the MNDWI, water bodies are maximally delineated from other classes. The EBBI, NDBI, and UI are better able to distinguish between built-up land and bare land, especially the EBBI. However, built-up land is confused with alpine meadows because they have similar values in the EBBI and NDBI. All bare land indices (NBLI, NDBAI, and NDSI) are able to enhance the bare land features. However, the NBLI enhances built-up land more than bare land, and the other landscape types are too close together in this index to be separated from each other. In the NDBAI and NDSI, landscape types overlap each other and cannot be separated except for water bodies and bare land.

Overall, water bodies are farthest from the rest of the landscape in the MNDWI, and bare land is farthest from the rest of the landscape in the EBBI. Due to the close proximity or even overlap of landscape response values to the NDBI, UI, NBLI, NDBAI, and NDSI, landscape types probably mix with each other in these indices and are not easily distinguished. Although bare land and alpine meadow can be confused in terms of the RVI and NDVI, they can be effectively distinguished by the WI and SVVI. The enhancement effect of the NDMI on water bodies is not as good as that of the MNDWI, but the NDMI has different response values for different landscape types due to moisture differences between landscapes. Finally, we selected seven feature indices for landscape classification: RVI, NDVI, WI, SVVI, MNDWI, NDMI, and EBBI.

2.3.4. Masking Out of Water Bodies and Selection of Stable Samples

In the MNDWI, the distance between the water body and the remaining landscapes becomes far enough to be easily separated from all others using the MNDWI. Based on the MNDWI, water bodies can be effectively separated from other landscapes (Figure 5), and the accuracy of water body extraction can reach more than 95% [31]. Thus, water bodies were extracted using a threshold value (0.6), and the mask of the water bodies was generated for each period. Then, the water bodies were masked out in the MVC of the RVI, NDVI WI, SVVI, NDMI, and EBBI images. These six feature indices and the DEM data (masked out water body area) were stacked together for each period.

Since there are little or no historical ground truth or other reference data, yet frequent cloud cover, it is difficult to obtain training (or validation) samples for each study period in this study region. To address these issues and to obtain consistent training (or validation) samples, a special method was used to derive these data from stable regions, places where the land cover does not change in the time series. Time-series SVVI images were employed to extract these regions, because the SVVI index can effectively distinguish various landscapes with water bodies removed (Figure 6). We used the linear trend (slope) of the MVC SVVI image over the time series to select the stability region. For any pixel x in the MVC SVVI image, the slope of that pixel on the time series is defined by Equation (1), and the Slope_svvi is calculated. The smaller the absolute value of the slope is, the less obvious the landscape change is. Therefore, we assume that pixel x is defined as a stable pixel if the Slope_svvi absolute value of the pixel is less than 0.05. The stable places with unchanged landscapes can be extracted in the time series [12].

(1)$Slop{e}_{svvi}=\frac{n\times {{\displaystyle \sum }}_{i=1}^{n}i{X}_i-{{\displaystyle \sum }}_{i=1}^{n}i{{\displaystyle \sum }}_{i=1}^n{X}_i}{n\times {{\displaystyle \sum }}_{i=1}^n{i}_{}^2-{({{\displaystyle \sum }}_{i=1}^{n}i)}^2}$

where i is the serial number of each research period from T1 to T16 (i = 1,2,3… 16), $X_i$ is the SVVI pixel value of the ith period, and n is the total number of research periods (n = 16).

Training (or validation) samples were extracted from the stable places. We created a stratified random sample of pixels from all stable areas to select the training (or validation) dataset, yielding approximately 3000 pixels for six landscapes (except water bodies) and between 300 and 600 pixels for each landscape type (except water bodies). Due to the lack of high-resolution images of the study area before 2009 in Google Earth, high-resolution images from 2009–2017 at two-year intervals in Google Earth were used as a reference to validate the precision of the extracted samples. The precision was more than 95%, and incorrect samples were eliminated. Then, 70% of the samples of each type were randomly selected as training pixels and the remaining 30% were selected as validation pixels for subsequent accuracy evaluation.

2.3.5. Landscape Classification and Landscape Dynamics Analysis

• (1). Landscape classification and accuracy assessment

Support vector machine (SVM) was chosen for this study because it is an efficient and robust supervised classifier that has been successfully applied to landscape classification [38]. Using the training samples, SVM was trained and used to classify the index stacking data into forest, alpine meadow, shrub, farmland, built-up land and bare land. Then, the water bodies were combined with the classification results.

For annual or biennial LTS maps over the long term, validating the accuracy is challenging as historical reference data or imagery are rare. Thus, the evaluation of LTS maps usually does not follow the model of standard accuracy evaluation and some compromise solutions are needed. To validate the LTS maps produced in this study, a compromise approach was adopted: we used a point-based accuracy assessment and assessed the accuracy of LTS maps using 30% of the sample from the stable places [12]. Since the water bodies are extracted through the MNDWI threshold with high accuracy [31], they are not included in the accuracy evaluation of the SVM classification results.

• (2). Analysis of landscape dynamics

To detect the landscape change information, we calculated the area and net changes (gains and losses) for each landscape type over consecutive periods. Cross-tabulation matrixes were generated to analyze the different landscape trajectories. We also estimated annual changes for each landscape type using the computational approach proposed by Puyravaud (2003) [39] (Equation (2)).

(2)$R=(\frac{100}{t_2-t_1})\times \ln (\frac{A_2}{A_1})$

where $A_1$ and $A_2$ represent the area of each landscape type at the beginning $t_1$ and the end $t_2$ of each period. $R$ is the percentage for each period.

In addition to landscape composition, the changing characteristics of landscape configuration were also identified. Four landscape pattern indices were calculated at the class level using FRAGSTATS v4.2: mean patch size (MPS), largest patch index (LPI), total edge length (TE), and mean Euclidean nearest-neighbor (MENN) distance. These four pattern indices were chosen because they provide key information to characterize and analyze the spatial configuration of landscapes [40,41]. A brief description of the four landscape pattern indices is shown in Table 1.

3. Results

3.1. Landscape Classification

Using the proposed five-stage framework, the biennial interval LTS maps were generated in the EHW, and the resulting landscape maps are shown in Figure 7. The overall accuracy of landscape classification was between 81.75% and 88.18% (Table 2), with the lowest value for T3 and the highest value for T16. For the individual landscape types, alpine meadow and bare land showed low values for both producer accuracy and user accuracy. Farmland was the landscape type with the highest classification accuracy, with values from 93.24% to 99.32% (user’s accuracy) and 89.57% to 98.64% (producer’s accuracy). Built-up land was also categorized as having good performance (user’s accuracy: 80.15–99.24% and producer’s accuracy: 94.62–99.23%). The accuracy values of the forest and shrub classes were intermediate, with values of 80.6–89.63% (producer’s accuracy) for the forest cover and 76.83–91.70% (producer’s accuracy) for the shrub cover (high value in only T16).

3.2. Changes in Landscape Composition

Forest was the predominant landscape throughout the entire study period, covering approximately 39.19% and 41.68% of the total study area at its beginning and the end, respectively. The total area of forest fluctuated during these years, but the overall trend increased (Figure 8). Forest cover increased moderately over the whole period, from 1138.82 km² to 1211.22 km², with an annual growth rate of 0.26%. Bare land decreased continuously during the 1987–2017 period, with a 46.07% net loss and an annual decline rate of 2.06%. The intermediate built-up area increased from 92.90 km² to 146.62 km², and the farmland area increased from 12.08% to 16.62% of the total area over the whole period. Built-up land and farmland continued to expand from T1, and the annual growth rates of the built-up land and farmland for the whole period were 1.52% and 1.06%, respectively.

Shrub areas showed a fluctuating trend over time, with a net gain of 1.26% for the total period. Water bodies decreased slightly over time from 8.96% to 8.88%, with a net change of −0.03%. Alpine meadow also shows fluctuating trends, but in different magnitudes over the whole study period from 9.3% to 8.26% (270.32 km² to 239.98 km²), with an 11.22% net loss and an annual reduction rate of −0.4% for the whole period.

Forest and shrub showed a very obvious state of mutual growth and decline in the study period of T4, T5, T6, T12, T13, T14, T15, and T16, that is, an increase in forest and a decrease in shrub or vice versa. Forest and shrub areas showed oscillating changes, mainly driven by climate change, as the regrowth of shrubs is often accompanied by an arid climate with little rainfall [42]. The areas most strongly affected by climate were located in valleys and high-altitude ridges.

3.3. Landscape Transformation

The transformation of the landscape types was captured for each period; the landscape response processes are shown in Figure 9 and Figure 10. The net changes in each landscape class over each research period are displayed in Figure 9. In particular, bare land showed the highest loss, with a net change by 228.97 km², followed by alpine meadow (net decrease by 30.33 km²) for the whole period (T1–T16). Farmland had the highest net gain (131.88 km²) for the whole study period, followed by forest (net gain 72.40 km²) and built-up land (net gain 53.72 km²). Shrubs and water bodies had the slightest net change within the whole study period (T1–T16); shrubs gained 3.73 km² and water declined by 2.42 km², respectively.

Figure 10 shows the transition process between one landscape class and the rest of the landscape classes in each study period. Over the period of T1–T16, built-up land was transferred in and out, and a large amount of new built-up land was developed on what was previously bare land, farmland, forest, and shrub, with most converted from bare land (19.11 km²), followed by farmland (17.76 km²), and 53.72 km² of land was converted to built-up land. Farmland area was gained from forest and bare land, and its lost area was converted into built-up land (Figure 10). For all study periods, a similar conversion phenomenon was observed, in which much bare land and forest were continuously converted to farmland. In general, farmland mainly increased. Forest gains occurred during T2–T3, T3–T4, T5–T6, T7–T8, T9–T10, T10–T11, T12–T13 and T14–T15, when the area was transformed to shrubland, alpine meadow, and bare land; the largest forest loss occurred during T13–T14, with a loss of 167.33 km², when the area was converted to mainly farmland, shrubland, and bare land (Figure 10). Throughout the study period, forest gained 142.63 km² and lost 70.24 km², for a net increase of 72.40 km². Alpine meadow mainly lost area to forest and gained area from two main landscape classes, shrub and bare land (Figure 10), over the entire study period. The greatest net loss of alpine meadow occurred in T13-T14, when it was largely converted into bare land (bare land gained 29.36 km²). The shrubland gained area from bare land and lost area to built-up land, farmland, forest, and alpine meadow (Figure 10). In the study period for T4–T5 and T5–T6, shrubs had a large area of net increase (137.69 km²) and net loss (168.58 km²). Bare land showed a loss (Figure 10). However, in T6–T7, T8–T9, and T13–T14, bare land was converted from built-up land, farmland, forest, and alpine meadow. The area of transition between water bodies and other landscape classes was very small (less than 5.50 km²) in all the study periods (Figure 10).

3.4. Changes in Landscape Pattern Indices

The landscape pattern differed significantly among the seven landscape classes (Figure 11). The high values of MPS and MENN for water were caused by the presence of more connected and compact water body patches (Figure 7). The LPI and TE index for water bodies remained almost constant, with the exception of MPS and MENN, which slightly decreased. The forest MPS index increased slightly in fluctuation, and its MENN index remained almost constant. After T13, forest showed a rapid increase in the LPI, while the LPI dropped to 10.82 at T16. Forests had higher LPI and TE values, indicating that forest patches had a larger area and more complex edge shapes. Built-up land showed a constant increase in the MPS index from 0.03 to 0.06. After T10 (LPI of 0.76), the LPI of built-up land started to increase rapidly, reaching a maximum of 2.00 at T15 and then decreasing to 1.80 at T16. The TE index of built-up land first increased and then decreased. In the study periods T7 and T11, the built-up land LPI reached a peak during two growth periods, which were 3209.31 and 3292.32, respectively. The MENN index for built-up land showed a moderate change, ranging from 161.47 to 185.19.

The LPI of farmland increased strongly from T13 to T16 and remained relatively constant in other periods; the highest LPI value was 11.11 at T14. The TE length of farmland increased from 3554.43 to 5516.76 with slight fluctuations. The MPS of farmland varied by less than 0.07 throughout the study period, with a maximum value of 0.19 at T3 and 0.12 at T9. The MENN of the farmland did not change drastically, and there was a very weak decreasing trend in the fluctuating changes. Alpine meadow showed a declining trend in the MPS between 1987 and 2017, with a change from 0.095 to 0.069. The alpine meadow LPI fluctuated significantly, reaching 2.72 at T6, decreasing to 1.07 at T7, and increasing to 2.01 at T9. The TE length of the alpine meadow exhibited a trend of first increasing and then decreasing, with a maximum of 4834.20 at T10. The MENN of the meadow also showed a trend of growth and then decline, reaching the maximum value (144.58) at T8.

The MPS of the shrubland consistently decreased, from 0.1 to 0.5, while the LPI increased in a fluctuating manner (from 0.80 to 3.12). The shrub TE length showed strong volatility, reaching a high value of 7311.57 at T9 and a low value of 2694.96 at T15. The MENN also showed obvious fluctuation characteristics, with a low value of 127.21 at T9 and a high value of 208.21 at T15. Notably, the TE length and MENN changes in shrubs showed opposite trends. Throughout the study period, landscape indices for bare land continued to decrease, with the exception of TE length, which had obvious fluctuations.

4. Discussion

Precise and long-term landscape dynamics are highly demanded for many applications of the EHW. However, the study area is located in a special geographical environment with cloudy coverage and variable topography, and there is a lack of historical reference data for obtaining classified samples. To address these challenges in extracting landscape information, we proposed an approach for obtaining landscape distribution data. We utilized this method to analyze the complex spatial and temporal changes in the landscape of the EHW. The landscape classification overall accuracy ranged from 81.75% to 88.18%, which met the accuracy requirements of most monographic studies. This result implied that the proposed method was effective and reliable in generating LTS maps for the EHW.

Based on the landscape classification results, forest was the predominant landscape throughout the analysis period, and forests dominated the eastern and western areas of the EHW. Alpine meadow also showed fluctuating trends. Shrub areas increased by 1.26% for the total period, and water bodies showed a small decrease with an overall net change of −0.03%. Built-up land and farmland continued to expand. Although the extent of increase in built-up areas was small (increased 53.72 km²), built-up areas can disproportionately influence the environment compared with changes in other landscape types. Bare land showed the highest loss, with a 228.97 km² net change. In the EHW, human activities and climate change have resulted in a gradual change in the landscape in the form of gains and losses. A comparison of the landscape classification and transformation analysis results indicated an increase in built-up land, farmland, and forest coupled with an obvious decrease in bare land over the entire study period. In the EHW, all landscape classes have changed or been converted into another class, and a substantially decreasing amount of bare land has taken place in the landscape.

The results showed a highly dynamic process for the landscape of the EHW; the region was under a continuous conversion between landscape types over the sixteen temporal periods. Based on LTS maps at biennial intervals, our study has illustrated the complex and varied landscape change processes that can occur in the EHW over a long period of time. Magnitude and trajectory changes in different landscape types occur over time, and such changes would hardly be detected by analyzing data obtained over coarse temporal intervals (e.g., 5–10 years). The high temporal resolution allowed us to understand in detail the linkage between the different change processes that occurred in the landscape; we want to highlight the importance of expanding the temporal resolution in landscape dynamics studies to effectively support sustainable environmental management.

5. Conclusions

To obtain information on landscape dynamics and to address the key issues encountered in extracting landscape maps (such as cloudy cover, lack of historical reference data), this study proposed a composite approach to generate biennial LTS maps and analyze landscape dynamics using Landsat time-series data in the EHW. We believe that the dense Landsat stacks methodology provides a successful pathway for quantifying landscape dynamics in a cloudy mountainous region by effectively minimizing the impacts of massive cloud cover [43]. Therefore, this study highlights the key role of Landsat data in the sustainable observation of Earth’s surface over a long time span. In the framework of our developed methodology, it should be emphasized that filling in missing data (cloud coverage, etc.) using the feature index MVC can help generate reliable time-series images for landscape extraction. Using of consistent stable samples can overcome the image difference problem. In particular, a set of stable samples used to train the classifier and validate the classification results was applied repeatedly for multiple periods without the need for reformulation [12,44]. This methodology proposed in this study is an important reference for similar research to be executed in other regions.

This study found that the landscape types in the EHW are dominated by forest, and the area of built-up land and farmland has expanded significantly; the area of water bodies did not show significant changes, and the area of bare land showed the highest loss throughout the study period. The landscape of the EHW has changed significantly, and the landscape dynamics are influenced by both regional social and economic development and climate change.

The analysis from the study revealed a detailed and complex landscape change process in the EHW. The findings of this study are significant in understanding the landscape dynamics in the cloudy mountainous regions, and provide scientific information for resource management, environmental protection, ecosystem service assessments, landscape planning and sustainability [45]. Hence, the study can help policy planners and makers to implement effective conservation and sustainable landscape management strategies in the EHW.

Author Contributions

Conceptualization, Y.C.; methodology, Y.C. and Y.Z.; software, Y.C.; validation, Y.C. and X.H.; formal analysis, Y.C., X.H. and J.F.; resources, Y.C. and Y.Z.; data curation, Y.C.; writing—original draft preparation, Y.C. and X.H.; writing—review and editing, Y.C., X.H., Y.Z. and J.F.; visualization, Y.C.; project administration, Y.C.; funding acquisition, Y.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Landsat imagery and DEM are openly available via the USGS.

Acknowledgments

We are grateful to the anonymous reviewers whose constructive suggestions have improved the quality of this study. We are grateful to the Yunnan Academician Workstation of Jingxiu Wang and the Plant Ecology Innovation Team of Dali University for their help. We are grateful for the data support from the ESPA (https://espa.cr.usgs.gov (accessed on 1 July 2019)). We also wish to express our gratitude to the Google Earth Pro platform for supplying the high spatial resolution maps for validating the accuracy of the stable samples.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Location and geographic extent of the study area.

Figure 2. The cloud-cover statistics for the study area from 1986 to 2017.

Figure 3. Temporal distribution (day of year) of Landsat images used in this study.

Figure 4. Flowchart of the proposed framework with dense Landsat stacks in the cloudy mountainous area.

Figure 5. The response of typical landscapes to the 12 indices in the study area. Here, the marks indicate the mean value of a landscape class.

Figure 6. SVVI for different landscape types (water bodies were removed).

Figure 7. Landscape maps for the start and end time periods, and detailed examples of landscape conversion over time.

Figure 8. Area change in different landscape types in the EHW from T1 to T16.

Figure 9. Net changes for each landscape type over each period.

Figure 10. Transition process between landscapes in the different periods.

Figure 11. Changes in the landscape pattern indices used for different landscape types in the EHW.

Table 1

Descriptions of selected landscape pattern indices.

Name	Unit	Description
MPS	km2	MPS is the average size of patch of a landscape type [41]. An increase in the MPS usually indicates a decrease in fragmentation.
LPI	%	LPI is a simple dominance measure, which is the percentage of the total landscape area constituting the largest patch [41].
TE	km	TE is the total edge length of a particular patch type [41]. Increasing the TE typically indicates an increase in the complexity of the patch shape.
MENN	m	MENN distance is the average shortest straight-line distance between nearest neighbor of the same class [41]. MENN has been used extensively to quantify patch isolation.

Table 2

Performance indicators for each landscape type (Pr: producer’s accuracy, Us: user’s accuracy).

Time	Overall Accuracy	Forest		Alpine Meadow		Shrub		Farmland		Built-Up Land		Bare Land
		Pr	Us	Pr	Us	Pr	Us	Pr	Us	Pr	Us	Pr	Us
T1: 1986–1987	82.72	81.87	85.49	76.66	78.57	82.00	79.46	92.00	93.24	98.20	83.21	74.78	78.18
T2: 1988–1989	81.99	83.01	80.13	78.08	81.43	76.83	77.13	95.33	96.62	97.22	80.15	70.54	82.73
T3: 1990–1991	81.75	81.35	83.91	72.82	77.50	85.15	75.58	89.57	98.65	96.52	84.73	73.21	74.55
T4: 1992–1993	82.80	81.82	82.33	73.38	76.79	83.75	77.91	96.67	97.97	99.19	93.13	72.27	78.18
T5: 1994–1995	82.40	81.42	82.97	77.27	78.93	78.63	75.58	92.31	97.3	100	93.13	73.39	72.73
T6: 1996–1997	83.12	83.97	82.65	72.35	80.36	84.28	74.81	96.05	98.65	98.39	93.13	74.14	78.18
T7: 1998–1999	82.32	80.6	85.17	73.00	78.21	83.7	73.64	95.92	95.27	96.77	91.60	75.68	76.36
T8: 2000–2001	86.17	84.33	84.86	80.63	81.79	81.99	82.95	96.67	97.97	99.23	98.47	86.00	78.18
T9: 2002–2003	85.29	86.08	83.91	78.05	80.00	79.70	82.17	96.00	97.30	99.22	96.95	84.62	80.00
T10: 2004–2005	84.00	81.87	85.49	73.67	73.93	82.79	78.29	98.64	97.97	99.23	98.47	81.98	82.73
T11: 2006–2007	83.92	83.44	82.65	74.59	80.71	82.55	75.19	98.00	99.32	99.22	97.71	76.99	79.09
T12: 2008–2009	84.08	87.50	81.70	78.28	81.07	78.36	81.40	89.94	96.62	95.97	90.84	82.24	80.00
T13: 2010–2011	87.30	89.63	84.54	78.38	82.86	87.11	86.43	95.97	96.62	98.48	99.24	80.36	81.82
T14: 2012–2013	85.61	87.5	83.91	78.69	81.79	83.97	85.27	95.3	95.95	94.62	93.89	78.70	77.27
T15: 2014–2015	86.98	87.03	86.75	77.74	83.57	88.46	80.23	97.33	98.65	99.23	98.47	80.53	82.73
T16: 2016–2017	88.18	85.85	88.01	80.34	83.21	91.70	85.66	96.05	98.65	97.73	98.47	85.58	80.91

References

1. Wang, Y.; Xia, T.; Shataer, R.; Zhang, S.; Li, Z. Analysis of Characteristics and Driving Factors of Land-Use Changes in the Tarim River Basin from 1990 to 2018. Sustainability; 2021; 13, 10263. [DOI: https://dx.doi.org/10.3390/su131810263]

2. Parihar, S.M.; Pandey, V.K.; Anshu,; Shree, K.; Moin, K.; Ali, M.B.; Narasimhan, K.; Rai, J.; Kamil, A. Land Use Dynamics and Impact on Regional Climate Post-Tehri Dam in the Bhilangana Basin, Garhwal Himalaya. Sustainability; 2022; 14, 10221. [DOI: https://dx.doi.org/10.3390/su141610221]

3. Yang, L.; Wang, Z.; Jin, G.; Chen, D.; Wang, Z. Geological risk assessment for the rapid development area of the Erhai Basin. Phys. Chem. Earth Parts A/B/C; 2015; 89–90, pp. 79-90. [DOI: https://dx.doi.org/10.1016/j.pce.2015.10.012]

4. Peng, J.; Ma, J.; Du, Y.; Zhang, L.; Hu, X. Ecological suitability evaluation for mountainous area development based on conceptual model of landscape structure, function, and dynamics. Ecol. Indic.; 2016; 61, pp. 500-511. [DOI: https://dx.doi.org/10.1016/j.ecolind.2015.10.002]

5. Peng, J.; Du, Y.; Liu, Y.; Hu, X. How to assess urban development potential in mountain areas? An approach of ecological carrying capacity in the view of coupled human and natural systems. Ecol. Indic.; 2016; 60, pp. 1017-1030. [DOI: https://dx.doi.org/10.1016/j.ecolind.2015.09.008]

6. Xu, J.; Wilkes, A. Biodiversity impact analysis in northwest Yunnan, southwest China. Biodivers. Conserv.; 2004; 13, pp. 959-983. [DOI: https://dx.doi.org/10.1023/B:BIOC.0000014464.80847.02]

7. Zhang, Y.; Wang, T.; Cai, C.; Li, C.; Liu, Y.; Bao, Y.; Guan, W. Landscape pattern and transition under natural and anthropogenic disturbance in an arid region of northwestern China. Int. J. Appl. Earth Obs. Geoinf.; 2016; 44, pp. 1-10. [DOI: https://dx.doi.org/10.1016/j.jag.2015.06.013]

8. Sulla-Menashe, D.; Friedl, M.A.; Woodcock, C.E. Sources of bias and variability in long-term Landsat time series over Canadian boreal forests. Remote Sens. Environ.; 2016; 177, pp. 206-219. [DOI: https://dx.doi.org/10.1016/j.rse.2016.02.041]

9. Li, X.; Gong, P.; Liang, L. A 30-year (1984–2013) record of annual urban dynamics of Beijing City derived from Landsat data. Remote Sens. Environ.; 2015; 166, pp. 78-90. [DOI: https://dx.doi.org/10.1016/j.rse.2015.06.007]

10. Nyland, K.E.; Gunn, G.E.; Shiklomanov, N.I.; Engstrom, R.N.; Streletskiy, D.A. Land Cover Change in the Lower Yenisei River Using Dense Stacking of Landsat Imagery in Google Earth Engine. Remote Sens.; 2018; 10, 1226. [DOI: https://dx.doi.org/10.3390/rs10081226]

11. Vuolo, F.; Ng, W.-T.; Atzberger, C. Smoothing and gap-filling of high resolution multi-spectral time series: Example of Landsat data. Int. J. Appl. Earth Obs. Geoinf.; 2017; 57, pp. 202-213. [DOI: https://dx.doi.org/10.1016/j.jag.2016.12.012]

12. Zhang, L.; Weng, Q. Annual dynamics of impervious surface in the Pearl River Delta, China, from 1988 to 2013, using time series Landsat imagery. ISPRS J. Photogramm. Remote Sens.; 2016; 113, pp. 86-96. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2016.01.003]

13. Moreira, E.P.; Valeriano, M.M. Application and evaluation of topographic correction methods to improve land cover mapping using object-based classification. Int. J. Appl. Earth Obs. Geoinf.; 2014; 32, pp. 208-217. [DOI: https://dx.doi.org/10.1016/j.jag.2014.04.006]

14. Li, H.; Wang, C.; Zhong, C.; Zhang, Z.; Liu, Q. Mapping Typical Urban LULC from Landsat Imagery without Training Samples or Self-Defined Parameters. Remote Sens.; 2017; 9, 700. [DOI: https://dx.doi.org/10.3390/rs9070700]

15. Griffiths, P.; van der Linden, S.; Kuemmerle, T.; Hostert, P. A Pixel-Based Landsat Compositing Algorithm for Large Area Land Cover Mapping. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2013; 6, pp. 2088-2101. [DOI: https://dx.doi.org/10.1109/JSTARS.2012.2228167]

16. Hermosilla, T.; Wulder, M.A.; White, J.C.; Coops, N.C.; Hobart, G.W. An integrated Landsat time series protocol for change detection and generation of annual gap-free surface reflectance composites. Remote Sens. Environ.; 2015; 158, pp. 220-234. [DOI: https://dx.doi.org/10.1016/j.rse.2014.11.005]

17. Xu, L.; Li, B.; Yuan, Y.; Gao, X.; Zhang, T. A Temporal-Spatial Iteration Method to Reconstruct NDVI Time Series Datasets. Remote Sens.; 2015; 7, pp. 8906-8924. [DOI: https://dx.doi.org/10.3390/rs70708906]

18. Malambo, L.; Heatwole, C.D. A Multitemporal Profile-Based Interpolation Method for Gap Filling Nonstationary Data. IEEE Trans. Geosci. Remote Sens.; 2016; 54, pp. 252-261. [DOI: https://dx.doi.org/10.1109/TGRS.2015.2453955]

19. Eivazi, A.; Kolesnikov, A.; Junttila, V.; Kauranne, T. Variance-preserving mosaicing of multiple satellite images for forest parameter estimation: Radiometric normalization. ISPRS J. Photogramm. Remote Sens.; 2015; 105, pp. 120-127. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2015.03.007]

20. Jiang, Z.-G.; Brosse, S.; Jiang, X.-M.; Zhang, E. Measuring ecosystem degradation through half a century of fish species introductions and extirpations in a large isolated lake. Ecol. Indic.; 2015; 58, pp. 104-112. [DOI: https://dx.doi.org/10.1016/j.ecolind.2015.05.040]

21. Wang, X.; Yu, S.; Huang, G.H. Land allocation based on integrated GIS-optimization modeling at a watershed level. Landsc. Urban Plan.; 2004; 66, pp. 61-74. [DOI: https://dx.doi.org/10.1016/S0169-2046(03)00095-1]

22. Guo, H.C.; Liu, L.; Huang, G.H.; Fuller, G.A.; Zou, R.; Yin, Y.Y. A system dynamics approach for regional environmental planning and management: A study for the Lake Erhai Basin. J. Environ. Manag.; 2001; 61, pp. 93-111. [DOI: https://dx.doi.org/10.1006/jema.2000.0400] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11381461]

23. Griffiths, P.; Kuemmerle, T.; Baumann, M.; Radeloff, V.C.; Abrudan, I.V.; Lieskovsky, J.; Munteanu, C.; Ostapowicz, K.; Hostert, P. Forest disturbances, forest recovery, and changes in forest types across the Carpathian ecoregion from 1985 to 2010 based on Landsat image composites. Remote Sens. Environ.; 2014; 151, pp. 72-88. [DOI: https://dx.doi.org/10.1016/j.rse.2013.04.022]

24. Nitze, I.; Grosse, G. Detection of landscape dynamics in the Arctic Lena Delta with temporally dense Landsat time-series stacks. Remote Sens. Environ.; 2016; 181, pp. 27-41. [DOI: https://dx.doi.org/10.1016/j.rse.2016.03.038]

25. Soenen, S.A.; Peddle, D.R.; Coburn, C.A. SCS+C: A modified Sun-canopy-sensor topographic correction in forested terrain. IEEE Trans. Geosci. Remote Sens.; 2005; 43, pp. 2148-2159. [DOI: https://dx.doi.org/10.1109/TGRS.2005.852480]

26. Qiu, S.; He, B.; Zhu, Z.; Liao, Z.; Quan, X. Improving Fmask cloud and cloud shadow detection in mountainous area for Landsats 4–8 images. Remote Sens. Environ.; 2017; 199, pp. 107-119. [DOI: https://dx.doi.org/10.1016/j.rse.2017.07.002]

27. Person, R. Remote mapping of standing crop biomass for estimation of the productivity of the short-grass Prairie. Proceedings of the Eighth International Symposium on Remote Sensing of Environment, Ann Arbor, MI, USA, 2–6 October 1972; Willow Run Laboratories: Brighton, MI, USA, 1972.

28. Tucker, C.J. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ.; 1979; 8, pp. 127-150. [DOI: https://dx.doi.org/10.1016/0034-4257(79)90013-0]

29. Lehmann, E.A.; Wallace, J.F.; Caccetta, P.A.; Furby, S.L.; Zdunic, K. Forest cover trends from time series Landsat data for the Australian continent. Int. J. Appl. Earth Obs. Geoinf.; 2013; 21, pp. 453-462. [DOI: https://dx.doi.org/10.1016/j.jag.2012.06.005]

30. Coulter, L.L.; Stow, D.A.; Tsai, Y.-H.; Ibanez, N.; Shih, H.-C.; Kerr, A.; Benza, M.; Weeks, J.R.; Mensah, F. Classification and assessment of land cover and land use change in southern Ghana using dense stacks of Landsat 7 ETM+ imagery. Remote Sens. Environ.; 2016; 184, pp. 396-409. [DOI: https://dx.doi.org/10.1016/j.rse.2016.07.016]

31. Xu, H. Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. Int. J. Remote Sens.; 2007; 27, pp. 3025-3033. [DOI: https://dx.doi.org/10.1080/01431160600589179]

32. Wilson, E.H.; Sader, S.A. Detection of forest harvest type using multiple dates of Landsat TM imagery. Remote Sens. Environ.; 2002; 80, pp. 385-396. [DOI: https://dx.doi.org/10.1016/S0034-4257(01)00318-2]

33. As-syakur, A.R.; Adnyana, I.W.S.; Arthana, I.W.; Nuarsa, I.W. Enhanced Built-Up and Bareness Index (EBBI) for Mapping Built-Up and Bare Land in an Urban Area. Remote Sens.; 2012; 4, pp. 2957-2970. [DOI: https://dx.doi.org/10.3390/rs4102957]

34. Zha, Y.; Gao, J.; Ni, S. Use of normalized difference built-up index in automatically mapping urban areas from TM imagery. Int. J. Remote Sens.; 2010; 24, pp. 583-594. [DOI: https://dx.doi.org/10.1080/01431160304987]

35. Hongmei, Z.; Xiaoling, C. Use of normalized difference bareness index in quickly mapping bare areas from TM/ETM+. Proceedings of the International Geoscience and Remote Sensing Symposium, Seoul, Korea, 25–29 July 2005; IEEE: New York, NY, USA, 2005.

36. Li, H.; Wang, C.; Zhong, C.; Su, A.; Xiong, C.; Wang, J.; Liu, J. Mapping Urban Bare Land Automatically from Landsat Imagery with a Simple Index. Remote Sens.; 2017; 9, 249. [DOI: https://dx.doi.org/10.3390/rs9030249]

37. Rogers, A.S.; Kearney, M.S. Reducing signature variability in unmixing coastal marsh Thematic Mapper scenes using spectral indices. Int. J. Remote Sens.; 2010; 25, pp. 2317-2335. [DOI: https://dx.doi.org/10.1080/01431160310001618103]

38. Zheng, B.; Myint, S.W.; Thenkabail, P.S.; Aggarwal, R.M. A support vector machine to identify irrigated crop types using time-series Landsat NDVI data. Int. J. Appl. Earth Obs. Geoinf.; 2015; 34, pp. 103-112. [DOI: https://dx.doi.org/10.1016/j.jag.2014.07.002]

39. Puyravaud, J.-P. Standardizing the calculation of the annual rate of deforestation. For. Ecol. Manag.; 2003; 177, pp. 593-596. [DOI: https://dx.doi.org/10.1016/s0378-1127(02)00335-3]

40. Walz, U. Indicators to monitor the structural diversity of landscapes. Ecol. Model.; 2015; 295, pp. 88-106. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2014.07.011]

41. McGarigal, K.; Marks, B.J. FRAGSTATS: Spatial Pattern Analysis Program for Quantifying Landscape Structure; U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station: Portland, OR, USA, 1995.

42. Cao, X.; Liu, Y.; Liu, Q.; Cui, X.; Chen, X.; Chen, J. Estimating the age and population structure of encroaching shrubs in arid/semiarid grasslands using high spatial resolution remote sensing imagery. Remote Sens. Environ.; 2018; 216, pp. 572-585. [DOI: https://dx.doi.org/10.1016/j.rse.2018.07.025]

43. Zhao, Y.; Feng, D.; Yu, L.; Cheng, Y.; Zhang, M.; Liu, X.; Xu, Y.; Fang, L.; Zhu, Z.; Gong, P. Long-Term Land Cover Dynamics (1986–2016) of Northeast China Derived from a Multi-Temporal Landsat Archive. Remote Sens.; 2019; 11, 599. [DOI: https://dx.doi.org/10.3390/rs11050599]

44. Stenzel, S.; Feilhauer, H.; Mack, B.; Metz, A.; Schmidtlein, S. Remote sensing of scattered Natura 2000 habitats using a one-class classifier. Int. J. Appl. Earth Obs. Geoinf.; 2014; 33, pp. 211-217. [DOI: https://dx.doi.org/10.1016/j.jag.2014.05.012]

45. Németh, G.; Lóczy, D.; Gyenizse, P. Long-Term Land Use and Landscape Pattern Changes in a Marshland of Hungary. Sustainability; 2021; 13, 12664. [DOI: https://dx.doi.org/10.3390/su132212664]