Evolution Simulation and Risk Analysis of Land Use Functions and Structures in Ecologically Fragile Watersheds

Abstract

The evolution of land use functions and structures in ecologically fragile watersheds have a direct impact on regional food security and sustainable ecological service supply. Previous studies that quantify and simulate land degradation in ecologically fragile areas from the perspective of long-term time series and the spatial structure of watersheds are rare. This paper takes the Huangshui Basin of the Qinghai-Tibet Plateau in China as a case study and proposes a long-time series evolution and scenario simulation method for land use function using the Google Earth Engine platform, which realizes the simulation of land use function and structure in ecologically fragile areas by space–time cube segmentation and integrated forest-based prediction. This allows the analysis of land degradation in terms of food security and ecological service degradation. The results show that: (1) the land use function and structure evolution of the Huangshui watershed from 1990 to 2020 have a significant temporospatial variation. In the midstream region, the construction land expanded 151.84% from 1990 to 2004, driven by urbanization and western development policy; in the middle and downstream region, the loss of farmland was nearly 12.68% from 1995 to 2005 due to the combined influence of the policy of returning farmland to forest and urban expansion. (2) By 2035, the construction land in the watershed will be further expanded by 28.47%, and the expansion intensity will be close to the threshold in the upstream and midstream areas and will continue to increase by 33.53% over 2020 in downstream areas. (3) The evolution of land use function and structure will further induce land degradation, causing a 15.30% loss of farmland and 114.20 km² of occupation of ecologically vulnerable areas, seriously threatening food security and ecological protection. Accordingly, this paper proposes policy suggestions to strengthen the spatial regulation for land degradation areas and the coordination of upstream, midstream, and downstream development.

1. Introduction

How human societies utilize, manage, and interact with land are key to addressing current sustainable development issues, including sustainable livelihoods, food security, biodiversity, climate change, and sustainable energy, which have been mentioned in high-level political agreements such as the 2030 agenda for sustainable development, the Paris climate agreement, and the convention on biological diversity [1,2]. Meanwhile, the global climate keeps showing a warming trend, which causes a long-term effect on the human–land system, especially in areas with vulnerable natural conditions [3,4]. Ecologically fragile areas are vulnerable to the combined effects of climate change and human activities due to the poor stability of their own functional and structural systems [5]. These areas are susceptible to exceeding critical thresholds and path dependence, leading to negative social or environmental impacts, and the impacts are difficult to reverse, are receiving more and more attention and becoming a hot topic. Retracing, monitoring, and analyzing the long-term time series changes of land use function and structure in ecologically fragile areas, revealing their process and mechanism, and modeling their potential degradation risks are of importance to global climate change monitoring and the sustainable development of ecologically fragile areas [6].

The existing dominant models for land use and land cover change include the logistic regression [7], Markov [8], system dynamics (SD) [9], cellular automata (CA) [10], CLUE-S [11], agent-based model (ABM) [12], and future land use simulation (FLUS) [13]. With the advantages of wide applicability, easy operation, and high openness to integration with other models, CA is widely used and has proven effective in various areas [14]. The key to these models is to set different development paths according to the land use evolution and development policies, predict the total amount and structure of land use in different paths, and use spatial simulation models to predict the land use patterns in the corresponding scenarios [13,15]. However, these methods are usually constrained by the temporal frequency of the historical land use data and lower classification accuracy, which lead to higher errors and uncertainties in prediction parameters, while seriously affecting the accuracy and reliability of simulation results.

With the open application of remote sensing big data platforms such as the Google Earth Engine, it is possible to use long-term free data and high-performance platforms for long-term land use classification [16]. At present, the realization of large-scale, long-term land use classification and change detection with the help of the Google Earth Engine (GEE) platform has become a hot issue, and a series of breakthroughs have been made. From a methodological view, these studies mainly include three categories: one is to identify and monitor the land use expansion through the construction of long-term remotely sensed ecological environment indexes such as NDVI [17], NDWI [18], etc. [19,20]; second is to realize the change detection of land use functions and structures [21,22] through full-element classification and mapping with long-term Landsat or Sentinel series images; third is to directly classify and extract change information by using multi-source images with the fusion method and dense stacking of time series data as input into the classifier, which prove that the accuracy of change detection can be significantly improved [23,24]. Based on long-term remote sensing classification results, some researchers explored the impacts caused by global issues, such as climate change [25,26], international protocols [27,28], and so on. However, these studies basically only focused on classification and change monitoring, and it is rarely shown how to reveal the evolution laws and patterns contained in long-term classification data and further apply them to the spatial simulation of future scenarios, which is our exploration in this study.

Located in the transition zone between the Qinghai-Tibet Plateau and Loess Plateau, the Huangshui Basin has fragile natural conditions and features of an alpine zone and soil erosion area, which causes great ecological vulnerability and significant ecological importance in the basin. Therefore, the Huangshui Basin is a typically ecologically fragile area. Meanwhile, it is the area with the highest population density and the most concentrated economy in the Qinghai-Tibet Plateau, and the human–land relationship is relatively adversarial and fragile. In 2020, the population of the Huangshui Basin was ~3.43 million, accounting for 56% of that of Qinghai Province; the population density was 95.62 persons/km², which was about 13 times that of Qinghai Province. In terms of agricultural production, the Huangshui Basin concentrates about 50% of it’s farmland in the Qinghai Province and produces about 65% of its grain; the industrial production there is more than 70% of that of the Qinghai Province. Under the background of global climate change, the changes in natural environment and the expansion of human activities in the past 30 years have had a relatively significant impact on the Huangshui Basin. In particular, the average annual precipitation shows a slight increase, the irrigation water diversion projects are continuously networked, the population and heavy chemical enterprises are further concentrated in the Huangshui Valley, the urban construction land has expanded by 3.6 times, and the tourist and economy in farming and pastoral areas of the upstream have continued to grow [29]. All these have led to great changes in the spatial structure of the Huangshui Basin, resulting in increasingly acute contradictions between urban construction, agricultural and animal husbandry production and ecologically functional land, and an imbalance in the spatial evolution of the upstream and downstream territorial structures. For example, the water pollution in the Xining-Minhe section was once in the lower five categories, and the downstream areas suffered from serious soil erosion and poverty problems [30]. However, China plans to achieve modernization by 2035, which means the intensity of human activities in Huangshui basin tends to maintain a high increase. In the future, with the development of the western region in the new period, the integration of Lanzhou-Xining urban agglomeration and the construction of the Qinghai-Tibet Plateau national park group, the population and economy of the Huangshui Basin will be further agglomerated, and the contradiction between humans and the land will become more prominent.

Therefore, this study took the Huangshui watershed of the Qinghai-Tibet plateau in China as an example and used the GEE platform to realize long-term land use classification. Specifically, a space-time cube that meets the prediction accuracy was constructed based on fully exploiting the advantage of long-term classification, the functional land parameters of future scenarios were measured and calculated by using space-time cube segmentation and prediction methods based on integrated forests, and the land-use function and structure in 2035 were simulated based on adaptive cellular automata. Based on the spatial simulation, relevant indicators were selected, and the risk of land degradation was analyzed from the two aspects of food security and ecological security. Finally, this study verified the reliability of the simulation and risk analysis results by combining field investigation visits to sit in typical areas of the watershed.

2. Study area and Data

2.1. Study Area

The Huangshui Basin (36°02′~37°28′N, 100°42′~103°04′E) is located in northeastern Qinghai, China, which is between the Qinghai-Tibet plateau and the Loess plateau. The basin comprises the mainstream area of Huangshui and the tributary Datong river basin, a plume-shaped water system with high terrain in the northwest and low terrain in the southeast. Its area is about 16,200 km² and its average elevation is 2300 m. The basin is divided into Chuanshui area, Qianshan area, Naoshan area and Shishan forest area from low to high, and the elevation of the dividing line is 2200 m, 2800 m, and 3200 m, respectively. From a sub-basin division perspective, the basin is divided into upstream, midstream and downstream and the points of demarcation are located in the confluences of Xina river and Yinsheng river (Figure 1), respectively. The climate here is a plateau arid and semi-arid continental climate, with high cold and little rainfall, and obvious vertical differences in climate and meteorology with an average temperature of 2.5~7.5 °C and large temperature differences between day and night. The average annual precipitation is 394.6 mm, about 70% of which is concentrated between June and September, resulting in frequent droughts and floods in the basin and a staggering trend. Water resources are scarce in the basin, the ecosystem is fragile, and droughts occur frequently, with typical characteristics of plateau basins and arid and semi-arid basins. The Huangshui basin involves Xining City, Haidong City, and Haibei Prefecture, including 13 districts.

2.2. Data

2.2.1. Remote Sensing Data

The remote sensing data applied in this study were mainly Landsat series image data, obtained from GEE collections “LANDSAT/LT05/C02/T1_L2” and “LANDSAT/LC08/C01/T1_SR”, including TM/ETM+/OLI, and periods are from June to October (vegetation flourishing period) and from November to March of the next year (vegetation withering period). In order to obtain as much observational data as possible with cloud cover <10% from 1987 to 2019, we utilized the pixel mosaic-based method [31] and >9000 images were obtained, obtaining 6 years images more than that of the traditional method. After image mosaic, image synthesis, and other operations, we obtained a long-term Landsat image dataset of the Huangshui watershed with a time step of 1 year. Spectral indices such as NDVI, NDWI, and NDBI were then calculated. Auxiliary data include VIIRS night light data (Nighttime Day/Night Band Composites Version 1) and digital elevation data (SRTM). The VIIRS nighttime light data was applied to distinguish urban and non-urban areas and is the monthly average radiation of nighttime data. The digital elevation data (SRTM) is obtained from the US Land Distributed Activity Archive and was gap-filled using open data (ASTER GDEM2, GMTED2010 and NED) [32].

2.2.2. Climate Data

The climate data in this study were FLDAS (Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System) data, including water content, humidity, evapotranspiration, soil temperature, and precipitation [33]. In addition, the Chinese meteorological dataset was obtained from 1915 stations in the Resource and Environmental Science and Data Center of the Chinese Academy of Sciences (http://www.resdc.cn/ (accessed on 7 July 2021)). After altitude correction, the annual average temperature (Ta), annual average precipitation (Pa), and accumulated temperature ≥ 10 °C (TaDEM), are obtained and then interpolated to obtain the nationwide dataset using the inverse distance-weighted average method.

2.2.3. Land Classification Sample Data

We selected a sample dataset of land-use function classification in the Huangshui basin year by year, including the data for cultivated land, shrubs, forests, grasslands, water bodies, urban construction, and unused land, by combining the visual interpretation of historical images and the selection of field samples. The selection of the dataset was mainly based on the 2019 high-resolution images of the same period as the selected Landsat image, combined with the field sample points to obtain the 2019 sample point data; after that, 2019 was used as the base year, and the reference was gradually selected forward. The image was interpreted sample point by sample point, and the sample points that changed the category were eliminated, and the corresponding sample points were filled to ensure the sample size. By repeating this process, a total of 1065 sample points covering 24-year were finally selected.

3. Methodology

The Methodology (Figure 2) in our study mainly included: (1) using a random forest and spatiotemporal consistency algorithm based on the GEE platform to realize long-term land use classification; (2) constructing a spatiotemporal cube and using spatiotemporal cube segmentation and integrated forest-based prediction methods to calculate the future scenario of the functional land parameters, and using the adaptive cellular automata to simulate land use function and structure; (3) selecting related indicators, such as the area and quality of food loss, the importance, and vulnerability of ecological service functions, analyzing the risks of land degradation from the perspective of food security and ecological security, and verifying the results with field investigations in typical areas of the watershed.

3.1. Long-Term Functional Classification of Land Use

In this study, we used the GEE platform and the random forest algorithm to initially classify long-term images year by year. We first obtained sample points of 2019 manually. Then to improve the accuracy and strengthen the consistency, we checked whether the features changed and provided category replacement based on historical images of 2018. The process was conducted annually until sample points of 1990 were obtained [29,34]. Next, the classifier was trained by continuously adjusting parameters such as sample point distribution, optimal feature vector combination, texture, optimal window size, and climate factors. After the training, we used the random forest classifier to classify the images of the corresponding years to obtain the initial classification results. The training was as follows: Firstly, taking 80% of the sample points as the training sample points and 20% as the verification sample points, and using the GEE’s built-in random function to randomly distribute the training sample points 30 times, we took the result with the highest accuracy as the final training sample [35]. Secondly, based on long-term remote sensing image data, the characteristic variables of each pixel, including the normalized vegetation index NDVI, normalized building index NDBI, normalized difference water index NDWI, and the night light data VIIRS, were calculated and combined year by year. Considering the significant influence of altitude and terrain factors on land-use function in the plateau area [26], the slope and aspect data generated by DEM were introduced as the classification basis, in which DEM adopted the SRTM dataset. At the same time, the parameters in the gray level co-occurrence matrix (GLCM) were used as texture feature variables, the window size and parameter combination with the highest verification accuracy were selected as the final texture feature, through the experiments of various window size and multi-type parameter combination. These texture feature variables included a total of 17 index parameters proposed by Haralick [36] and Conners [37], which are angular second moment (ASM), correlation (CORR), entropy (SENT), variance (VAR), inverse difference moment (IDM), average (SAVG), sum variance (SVAR), sum entropy (SENT), information entropy (ENT), difference variance (DVAR), difference entropy (DENT), relevant information measure 1 (IMCORR1), relevant information measure 2 (IMCORR2), maximum correlation coefficient (MAXCORR), dissimilarity (DISS), inertia (INERTIA), cluster shadow (SHADE), and cluster protrusion (PROM). In addition, considering the natural background characteristics of the Huangshui basin located on the Qinghai-Tibet plateau is sensitive to hydrothermal conditions, meteorological elements were introduced for classification. Among them, the temperature was characterized by annual average temperature, annual accumulated temperature and soil average temperature, and the precipitation was characterized by annual precipitation, soil moisture and soil moisture content. Finally, the final training samples, feature variables, and terrain, texture and climate features were input into the random forest classifier.

Due to the interannual differences in the quality of the Landsat images used in the long-term classification results, we implemented a two-way consistency detection algorithm by constructing a sliding window and a seed window to correct the consistency of the classification results [29]. The main processing included: (1) setting up the bidirectional sliding window $w_m$ and the seed window $w_s$ by selecting the five years with the highest average precision as the starting position of the sliding window, and the third year as the reference year of the seed window to construct a 3 × 3 seed window; (2) placing the sliding window $w_m$ immediately outside the seed window $w_s$; (3) determining the dominant classification $f_m$ in the sliding window $w_m$ by investing the frequency ρ of a function classification: if it was >0.7, the function classification was regarded as the dominant classification $f_m$ of the sliding window $w_m$, and if the classification frequency of all functions was <0.7, it was considered that there was no dominant classification in the sliding window; (4) determining whether the object type fs in the seed window $w_s$ was consistent with the dominant classification $f_m$ in the adjacent sliding detection window $w_m$: if they were consistent, we set all the classifications in $w_m$ to $f_s$, and moved the position of the seed window to the grid position where the last dominant classification $f_s$ appeared in the sliding window, and went to step (2) until the seed window reached the beginning and end of the time series; if it was inconsistent, we moved the seed window one bit inward, and went to step (2), and repeated until the seed window reached the beginning and end of the time series to complete the consistency correction.

3.2. Scenario Simulation of Land Use Function

Simulation analysis of land use function mainly included: (1) land use scale simulation based on grid division and random forest, (2) spatial structure simulation based on adaptive cellular automata, which is the simulation of the spatial structure of land use functions in 2035 based on land use function classification data, quantitative impact factor parameters and the result of land use scale simulation.

The Huangshui basin had the characteristics of significant spatial differentiation of natural background conditions, the factors that dominated the evolution of land-use functions in different regions were significantly different, and it was difficult to obtain simulation results with high confidence. Therefore, we constructed a 2 km × 2 km grid division unit, downscaled the long-term land use function classification data, and calculated the area proportion of each functional type $P_t^{m,q}$ as a function parameter year by year. The $P_t^{m,q}$ was expressed as follows:

$$P_t^{m,q}=\frac{s_t^{m,q}}{s}$$

where $s_t^{m,q}$ represents the area occupied by the land-use function q in the grid unit m at the t-th time, and s is the area of the grid unit.

Then, based on the existing functional parameter data, the Lagrangian interpolation method was used to interpolate the functional parameters of 8 time nodes that lacked land use functional classification data grid by grid, and the long time series data of functional parameters with a step of 1 year were constructed. After that, based on the long-term time series data of grid cell functional parameters, the random forest algorithm was used to predict the functional parameters for 2035. We implemented the algorithm through the tool Forest-based Forecast, developed and provided in ArcGIS Pro 10.8 by ESRI. The size of the decision tree forest constructed in this study was 300, and the length of the decision tree training sample was 7 time nodes. To improve prediction accuracy, the data of three time nodes were selected as the verification samples, and the decision tree with the root mean square error (RMSE), calculated by simulated value and observed value, less than 10⁻⁴ was retained by simulating the verification samples based on the training samples. Finally, the simulation results of functional parameters were normalized by functional classification, and the total area of each functional classification in the Huangshui basin was counted to obtain the simulation results of land use functional classification in the Huangshui basin in 2035.

Based on the principle of functional genetics, we selected corresponding indicators by considering natural carrying capacity, development status, and future potential. Critical indicators such as land construction suitability, water resource convenience, and geological disaster risk were selected for natural carrying capacity. Critical indicators such as population density and economic agglomeration were selected for development status. The future potential focused on the quantification of location conditions, with traffic dominance as the main indicator including traffic accessibility and distance from the central city. Then, based on the attributes of functional classification, the quantitative parameters of different influencing factors were linear combination by category to calculate the functional suitability probability $S_k$ of land-use. The functional suitability probability $S_k$ of land-use was expressed as follows:

$$S_k={{\displaystyle \sum }}_x{w}_{x,k}\ast I_x$$

where x is the impact factor, $w_{x,k}$ is the impact weight of the impact factor x on the functional classification k, and $I_x$ is the quantitative parameter of the impact factor x.

Using $S_k$ as the driving factor, and the land scale simulation result of functional classification as the target parameter, the generation probability $T{P}_{p,k}^t$ of each functional classification was calculated by adaptive cellular automata, which is conducted in GeoSOS-FLUS V2.3, a software developed in 2017 [13], and the evolution simulation of functional classification was realized by generating floating-point random numbers [13].

$$T{P}_{p,k}^t=S_{p,k}\ast {\Omega }_{p,k}^t\ast Interti{a}_k^t$$

where k is the land type; p is the spatial location; t represents the t-th iteration; $T{P}_{p,k}^t$ represents the probability that the functional classification; $S_{p,k}$ is the functional suitability probability of land-use; ${\Omega }_{p,k}^t$ represents the neighborhood influence; and $Interti{a}_k^t$ represents the inertia coefficient which is used to realize the adaptive process of the cellular automaton.

3.3. Risk Analysis of Land Degradation

We analyzed the risk of land degradation from the aspects of food security and ecological security in this study. Food security risk mainly refers to the loss of farmland quantity and quality due to the superposition of human activities and climate change, which is mainly characterized by landscape indexes such as total farmland, farmland quality, average patch area, the proportion of large-scale patches, and fragmented patches in total patches. Ecological security risk refers to the decline of ecological service functions and the enhancement of ecological vulnerability due to human activities and climate change, which refers to the encroachment of other land-use on the ecological service function land, mainly characterized by the importance evaluation of ecological conservation. The importance evaluation (Z) was expressed by ecosystem service function importance (X) and ecological vulnerability (Y), both of which were measured by the quality index [38], InVEST model [39] and so on. Among them, the X was related to the importance of ecosystem services such as water conservation [40], soil and water conservation [41], biodiversity maintenance [42], and windbreak and sand fixation [43]; the Y took soil erosion [44] and desertification vulnerability [45] into consideration:

$$\left\{\begin{array}{l}Z=f\left(X,Y\right)\\ X=C_x\left(x_1,x_2,x_3,\dots ,x_n\right)\\ Y=C_y\left(y_1,y_2,\dots y_n\right)\end{array}\right.$$

where Z represents the evaluation level of ecological protection importance, X and Y represent the importance of ecosystem service functions and ecological vulnerability, respectively, $x_1$~$x_n$ represent different types of ecological service functions, and y₁~y_n represents different vulnerability types.

4. Results

4.1. Long-Term Land Use Functional Classification Result

The long-term classification results of land use in the Huangshui Basin are shown in Figure 3. The classification results show that the overall accuracy of the initial classification results was 76.02~88.92%, and the average accuracy was 83.06% with an improvement of 2.32%. Grassland and farmland were the main land-use types in this basin, which accounted for >77% of the total area of the basin. Grassland was mainly distributed in alpine mountain areas and river valley areas, among which the area of river valley grassland always remained higher than 4652.14 km². The farmland was mainly distributed in the river valley area, among which the area with a slope of <12° in the basin reached 2778.96 km² in 2019, accounting for 55.98% of the total area. The woodland was mainly (63.60%) distributed in mountainous areas with higher altitudes. Compared with grassland, woodland areas were generally located upstream with higher elevations and steeper slopes.

The spatial pattern of land use in the basin varied not only with terrain, but also with upper, middle, and lower streams. The land-use types dominated by human activities were mainly distributed in the middle and lower steams, while the original natural land types were mainly in upstream. In 2019, the farmland area in the midstream and downstream reached 3658.18 km², accounting for 73.70% of the total area of the watershed; for construction land dominated by human activities as well, it mainly (725.88 km², 62.72%) located in Xining and Haidong in midstream area. The bare land was mainly distributed upstream, with an area of 348.47 km², which was higher than the sum of that in the middle and lower stream, accounting for 57.18%.

We found that the largest changes were in construction land and farmland from 1990 to 2019. The construction land increased from 315.72 km² to 1091.83 km², with an expansion of 2.67 times; while the area of farmland decreased in a wide range, from 6084.17 km² to 4900.33 km², with a loss of 18.41%, which is the land type that lost the largest areas over the past 30 years. A major reason for the loss of farmland was the expansion of construction land, occupying a total area of 601.73 km², accounting for 77.53% of the total area of newly added urban construction land; another major reason for this was that the evolution to ecological land and the area of woodland, shrubs and other ecological land from evolution was 762.21 km², which was higher than the amount of urban construction encroaching. At the same time, the farmland showed a certain trend of migrating to the mountainous areas, which was revealed by the average slope of the farmland in 2019 as 12.08° and as 11.86° in 1990.

In addition, the land use change in this watershed had typical stage characteristics. First, the rate of urban expansion in different stages was significantly different and the main new urban construction land was concentrated in the period from 2000 to 2014, accounting for 56.54% of the total new area. Among them, from 2000 to 2010, the average annual growth rate was 3.84%; from 2010 to 2014, it was a high-speed expansion stage, with an average annual growth rate of 6.74%, which was more consistent with the region’s response to the western development strategy. Second, the process of encroaching on farmland by construction land and forests and grasses had two stages: 2000–2006 and 2014–2018, with 362.30 km² and 160.21 km² of farmland being converted into woodland respectively, which is closely related to the project of “returning farmland to forest”.

4.2. Scenario Simulation Result of Land Use Function

The predicted parameters of different land use functions and structures in the Huangshui basin for 2035 are shown in

Table 1. Simulation scale result table (km²) + specific gravity.

	Arable	Shrubs	Woodland	Grassland	Water	Construction	Bare Land
2019	4966.18	2231.49	568.51	10524.76	60.53	1157.88	609.72
2035	4206.25	2509.25	455.31	10611.83	96.66	1487.57	752.20

. During the period from 2020 to 2035, the construction land in this basin still maintains a trend of rapid expansion. The simulation results reached 1487.57 km², which increased by 330.31 km² compared to 2019, and the average annual expansion rate reached 1.69%. The ecologically degraded bare land is another type of area that increased significantly, which increased from 609.72 km² to 752.20 km², with an increase of 23.37%. Shrubs and water showed a slight expansion trend, expanding by 277.76 km² and 36.13 km², with an increase of 11.07% and 37.38%, respectively. However, a loss of farmland in the simulation is always serious. The simulated farmland was 4206.25 km², which was only 84.70% of that in 2019, and the average annual decline rate reached 1.10%. The woodland also had serious degeneration, which decreased from 568.51 km² in 2019 to 455.31 km², with a decline rate of 19.91%. The grassland did not change significantly, and the area change was only 87.07 km², accounting for 0.82% of the total area.

The land use spatial simulation of the Huangshui basin for 2035 with the adaptive cellular automata model is shown in Figure 4. From the perspective of terrain structure, there were obvious differences in land-use in this Basin. The construction land and farmland are concentrated in river valleys and shallow mountainous areas with an altitude of <2800 m, with an area of 1280.01 km² and 2938.26 km² respectively, accounting for 86.20% and 69.40% of their respective total area. The woodlands were mainly distributed in the origins of the main stream and main tributaries of the Huangshui river and the upper mountains, and were concentrated in the water conservation areas with an altitude of >2800 m. Shrubs were mainly distributed in the outer edge of the valley, concentrated on the edge of farmland; grassland was widely distributed in the whole basin; and bare land was densely distributed in the mountainous areas of the watershed or part of the valley slopes of the tributaries of Huangshui.

From the perspective of the main and tributary network of the river, the simulation results had obvious spatial differences between the upstream and downstream and the main and tributary. The overall function and structure of the upper streams of the watershed had a relatively small evolution, and were dominated by natural land such as forests and grasses. At the same time, there was also the expansion of urban construction land and the loss of farmland, but the evolution scale was small. The area of newly added urban construction land and the area of lost farmland were 61.38 km² and 121.33 km², with ranges of 23.24% and 9.29%, which were lower than the overall level of the basin. The middle stream had the areas with the highest intensity of human activities and growth rate, and construction land and farmland accounted for 51.75% and 56.64%. From 2019 to 2035, the new construction land area was 167.93 km², accounting for 50.85% of the total new expansion; in addition, the mountain bare land at the source of the tributaries in the middle stream had the most obvious expansion, with an increase of 81.19 km², accounting for 50.69% of the new bare land area, and the expansion rate reached 163.36%. The expansion of construction land in the lower stream is also extremely serious. The newly added urban construction land covered an area of 100.88 km², with an average annual growth rate of 1.95%, which was 0.25% higher than the average level. Similar to the midstream, the bare land in downstream had a substantial expansion as well, with an expansion area of 38.05 km² and an expansion rate of 54.54%.

4.3. Risk Analysis Results of Land Degradation

The simulation results of land use in the Huangshui basin demonstrated that both construction land and bare land showed a significant expansion trend, which was accompanied by the loss of a large number of high-quality farmland and the loss of ecological service functions, threatening the food security and biodiversity of the watershed. The area of farmland decreased significantly, and the degree of fragmentation continued to increase. By 2035, 14.70% of the farmland was lost in this basin, and the average patch area was 0.326 km², with a decrease of 4.66%; and the broken farmland with an area of <1000 m² increased from 64.44% to 74.64%. By analyzing the evolution from farmland to other land use types (Figure 5), the expansion of construction land (typical of Duoba Town shown in Figure 5b), the conversion between woodland and grassland and farmland, and the degradation of farmland to bare land (typical of Hetaozhuang Town shown in Figure 5c) were the main reasons for the loss of farmland, accounting for 82% of the total farmland loss. The farmland replaced by woodlands and grasslands was mainly distributed in the sloping farmland with an altitude of >2800 m and a slope of >25% with a lower quality level. In addition, the degradation of farmland to bare land may also occur in the sloping areas of river valleys.

In addition to the loss of farmland, the risk of land degradation was also reflected in the partial weakening and loss of ecological services (Figure 6). The importance evaluation of the ecological service functions showed that the area with moderate or above ecological service function importance in the basin is 13,440.65 km², accounting for 66.78% of the total area of the basin. Among them, the proportion of soil and water conservation and biodiversity was relatively high, accounting for 88.09% and 97.37%, respectively, mainly distributed in the valley edge of this basin. This area is relatively flat, and the water and heat conditions are relatively abundant, which is suitable for the growth of various plants, and plays an important role in maintaining regional biodiversity and soil and water conservation. By 2035, the loss of areas with important ecological service functions increased to 111.64 km² (typical in Ledu District shown in Figure 6b), with an increase of 67.43% over the base year. The lost areas were mainly distributed in the transition zone from the middle and lower streams, accounting for 28.71% and 19.47% of the total area of important ecological service functions, respectively. The ecological vulnerability analysis showed that the total area above moderately vulnerable in this basin was 6792.41 km², accounting for 33.75%, mainly concentrated in the valley slopes and the watershed mountain areas of the basin. The area of bare land increased significantly, and the natural land degradation was relatively serious. Among them, 114.20 km² of newly added bare land was in areas with strong ecological vulnerability (typical of Lierbao Town shown in Figure 6c, where 52.99% bare land expansion was simulated in the area with strong ecological vulnerability), accounting for 71.29% of newly added bare land. The expansion of bare land in ecologically fragile areas mainly occurred on the slopes of tributary valleys in the middle stream. At the same time, by superimposing the results of the individual evaluation of ecological vulnerability and the simulated expansion of bare land, we found that 70.35% of the newly developed bare land was located in the middle and high-risk areas of soil erosion in this basin.

5. Discussion

The data used for the classification of land use functions in this study were long-series Landsat satellite images. Due to the aging of Landsat 5 sensors and the serious banding of Landsat 7 satellite images from 2007 to 2013, there are certain flaws in the reliability of the data [46]. The advantages of the method in this paper have been confirmed in previous studies, and a better classification has been achieved [47,48]. The overall accuracy (OA) of Random Forest classification results ranges from 76.38% to 89.07% and the Kappa range from 0.6881 to 0.9131. After bi-detection process, the OA averagely improved 5.60%, ranging from 84.20% to 91.74%.

Meanwhile, the method in this paper combines random forest (RF) and cellular automata (CA), both of which have proven stable and robust in land use research [49,50,51,52]. Moreover, a series of works have conducted a sensitivity analysis of the CA model, which found that the results are sensitive to time step, impact of driving factors, and neighborhood influence [13,53]. Because the CA model is not sensitive to the area that random forest simulates, the combined method is also stable and robust. To further assess accuracy, we used the indicator “figure of merit” (FoM) as a veridiction method. FoM is a ratio representing the proportion of samples that are observed and simulated correctly [54] and has been widely used in previous studies [55]. We simulated land use in 2018 using our method and the Markov-FLUS model, separately. The FoM using the method in this paper was 12.63%, which was 1.33% higher than FoM using Markov-FLUS model.

As a high-intensity human activity area with an important but vulnerable environment, the Huangshui basin is a hotspot of environmental protection policies and development plans, such as “returning farmland to forest”, “protection of prime farmland”, “Lanzhou-Xining Urban Area Development Plan”, and so on. Effected by multiple policies, interconversion of land use has occurred in the Huangshui Basin, between farmland and shrubland, for instance. Taking policy factors and multiple special types of areas, such as agricultural and animal husbandry interlaced zones and ecologically fragile areas into consideration, the driving force is relatively complex, and the reliability of prediction and simulation needs to be verified. Therefore, this study conducted a field survey, and verified the conclusions of this study through a comprehensive practical investigation and visual interpretation of historical image data for the typical sites involved in the study.

Because of the loss of cultivated land and ecological degradation caused by the expansion of construction land, we selected Duoba Town in Huangzhong District Figure 7a and the urban area of Ledu District Figure 7b as the verification and inspection sites, both of which are urban expansion areas in the “Lanzhou-Xining Urban Area Development Plan”. Combining historical images and land use function classification data, it was found that Duoba Town experienced a rapid expansion of urban construction land after 2012, and the loss of its cultivated land accelerated during the same period; in the simulation results, Duoba Town maintained the rapid expansion of urban construction land, its area reached 58.44 km² in 2035, with an increase of 15.73 km² compared to that in 2019; among them, 14.01 km² of new urban construction land was converted from cultivated land, accounting for 89.08% of the totally new area. The Ledu District had a similar situation to that of Duoba Town, and the overall urban area along the main banks of the Huangshui river showed a significant increase in 2013. The area of newly added urban construction land was 7.60 km², of which the area of the area with medium or above ecological service function importance accounted for 44.38%. The new urban area of Ledu District is flat, which is convenient for the construction implementation, and the transportation infrastructure is relatively superior, which reveals the extremely strong driving force for evolution into urban construction land. Considering that Ledu District is located at the entrance of the tributaries in the low stream, the weakening of the soil and water conservation service function can easily lead to the instability of the regional ecosystem, resulting in the risk of soil erosion in the upstream area.

Given the degradation of farmland and ecological area into bare land, we selected Hetaozhuang town in Minhe County Figure 7c and Donggou Town in Huzhu County Figure 7d as verification sites. Several protection projects were conducted on these sites, such as “returning farmland to forest” and “protection of prime farmland”. Based on the classification data of land use function, it was found that Hetaozhuang Town continued to maintain a certain decline of cultivated land after 1990, and the bare land also expanded slightly. In the simulation results, the cultivated land in Hetaozhuang Town continued to decline at significantly increased rates, which mainly evolved into grassland and bare land. The area of cultivated land in decline reached 4.91 km², accounting for 20.41% of the total cultivated land. From the field investigation, it was found that there is a significant vertical difference in some mountain land types in this area, which shows that the top mountain is land or bare land, and the bottom is cultivated land and vulnerability may induce bare land to encroach on cultivated areas. The bare land in Donggou Town also showed a large-scale expansion in the simulation results, with an expansion area of 12.24 km², and the simulated degradation area was 9.07 km², revealing potential risks of land degradation in areas with strong ecological vulnerability. The investigation shows that the vegetation coverage in this area is relatively good, the main natural vegetation is grassland, and the artificial development of the mountain is common with cultivated land at the bottom of the mountain. Further investigation demonstrated that this area consisted of bare land and low-coverage grassland before 2015, and there was an expansion trend, and the implementation of the wasteland treatment project and the conversion of farmland to forest has gradually restored the wasteland and part of the cultivated land into grassland and forest land after 2015. The stability of the ecosystem in Donggou Town has been enhanced, which has greatly curbed the expansion trend of bare land. With the auxiliary UAV images, we found that the grassland in Donggou town has regional heterogeneity and large-scale differences in coverage levels, which can prove that there are ecological restoration projects carried out in stages in this region.

After investigation and verification of typical sites in the Huangshui basin, it was found that: (1) there is a strong driving force for land-use conversion in each area based on the verification of the background conditions of urban construction land expansion and bare land degradation areas, revealing that the evolution of functional land is reasonable and the conclusions are reliable; (2) the degradation of bare land in the watershed mainly result from the regional topography, soil quality and vegetation coverage, the change of which has a great influence on the regional ecological effect. Destructive human activities lead to irreversible land degradation and bring serious ecological risks, while ecological management projects such as returning farmland to forests can curb the trend and improve regional ecological security.

6. Conclusions

This paper modeled the spatial evolution of land use function and structure in the Huangshui basin on the Qinghai-Tibet plateau using long-term time-series remote sensing land use/land cover classification data and analyzed land degradation analysis in terms of protection and degradation of ecological service functions. The study shows that Huangshui basin, as an ecologically fragile watershed, has undergone drastic changes in the past 30 years, manifested by the massive expansion of construction land and the massive loss of farmland. This change has been influenced not only by human activities such as urbanization and ecological immigration, but also by the combined effects of various policies such as the return of farmland to forests, farmland protection, and climate change. By continuing to follow past development paths, the Huangshui basin will suffer from further extensive land degradation, as evidenced by the massive expansion of built-up land and bare land, which will encroach further, leading to the loss of farmland, as well as the area of high ecological importance and vulnerability. Given the spatial management and control of specific functional areas and the coordination of upstream, midstream, and downstream development, we put forward specific policy recommendations in this study: (1) strengthen the protection of cultivated land, delineate the farmland protection redline, promote scientific planting pattern, and encourage the development of water-saving agriculture to improve planting efficiency; (2) strengthen ecological protection, ensure the ecological service function of the river basin, delineate ecological protection redline, ensure the ecological service function of the river basin, maintain the stability of the river basin ecosystem, and maximize the ecological benefits; (3) promote the unified planning and improve the collaborative management of the river basins.