**Accessing the Time-Series Two-Dimensional Displacements around a Reservoir Using Multi-Orbit SAR Datasets: A Case Study of Xiluodu Hydropower Station**

**Qi Chen 1,2, Heng Zhang 1,2, Bing Xu 3,\*, Zhe Liu <sup>4</sup> and Wenxiang Mao <sup>3</sup>**


**Abstract:** The construction of large-scale hydropower stations could solve the problem of China's power and energy shortages. However, the construction of hydropower stations requires reservoir water storage. Artificially raising the water level by several tens of meters or even hundreds of meters will undoubtedly change the hydrogeological conditions of an area, which will lead to surface deformation near the reservoir. In this paper, we first used SBAS-InSAR technology to monitor the surface deformation near the Xiluodu reservoir area for various data and analyzed the surface deformation of the Xiluodu reservoir area from 2014 to 2019. By using the 12 ALOS2 ascending data, the 100 Sentinel-1 ascending data, and the 97 Sentinel-1 descending data, the horizontal and vertical deformations of the Xiluodu reservoir area were obtained. We found that the Xiluodu reservoir area is mainly deformed along the vertical shore, with a maximum deformation rate of 250 mm/a, accompanied by vertical deformation, and the maximum deformation rate is 60 mm/a. Furthermore, by analyzing the relationship between the horizontal deformation sequence, the vertical deformation sequence, and the impoundment, we found the following: (1) Since the commencement of Xiluodu water storage, the vertical shore direction displacement has continued to increase, indicating that the deformation caused by the water storage is not due to the elastic displacement caused by the load, but by irreversible shaping displacement. According to its development trend, we speculate that the vertical shore direction displacement will continue to increase until it eventually stabilizes; (2) Vertical displacement increases rapidly in the initial stage of water storage; after two water-storage cycles, absolute settlement begins to slow down in the vertical direction, but its deformation still changes with the change in the storage period.

**Keywords:** SBAS; two-dimensional deformation sequence; water-level adjustment

#### **1. Introduction**

As a clean and recyclable resource, hydroelectricity is receiving increasing attention from all countries [1]. Inevitably, the rapid development of China's economy has led to a huge demand for energy. Southwest China has abundant water resources for the construction of hydropower stations, so it has obtained a reputation as the "Asian battery" [2]. However, the impounding of a reservoir changes the conditions of local hydrogeology and the regional gravity field, resulting in ground deformation [3]. Furthermore, fluctuations in water level accelerate the deformation, causing local geohazards. In 1961, a landslide with a volume of about 165 × <sup>10</sup><sup>4</sup> m3 occurred in the Zixi Reservoir in Zishui, Hunan Province, causing 40 deaths [4]. In 1963, a landslide occurred in the Vajont Reservoir in Italy; a huge volume of landslide soil (about 270 × 106 <sup>m</sup>3) rushed into the reservoir, and the surge caused 1925 deaths [5].

**Citation:** Chen, Q.; Zhang, H.; Xu, B.; Liu, Z.; Mao, W. Accessing the Time-Series Two-Dimensional Displacements around a Reservoir Using Multi-Orbit SAR Datasets: A Case Study of Xiluodu Hydropower Station. *Remote Sens.* **2023**, *15*, 168. https://doi.org/10.3390/rs15010168

Academic Editors: Qiusheng Wu, Xinyi Shen, Jun Li and Chengye Zhang

Received: 15 November 2022 Revised: 12 December 2022 Accepted: 19 December 2022 Published: 28 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The Three Gorges Dam is the largest hydropower station in the world. After the water level of the Three Gorges Dam reached 175 m in 2008, the number of resurrection landslides was tens of times the number from 1992 to 1995 [6]. In order to decrease such geological disasters, scholars are paying more attention to the ground deformation around hydropower stations.

Currently, the tools to monitor ground deformation include the spirit leveling method [7], the Global Positioning System (GPS) [8], and optical remote sensing [9]. However, the spirit leveling and GPS methods are costly and a have low spatial resolution for capturing the detail of a whole deformation pattern [10]. Optical remote sensing is susceptible to the air conditions of clouds and rain.

Synthetic aperture radar interferometry (InSAR), a satellite-based geodetic technique, can overcome the influence of cloud and rain on optical remote sensing and can work both day and night, demonstrating advantages in ground-deformation monitoring [11]. InSAR has been used in the geohazard monitoring of landslides [12], earthquakes [13], volcanos [14], and ground subsidence [15], among others. Additionally, in order to bypass the effects of decorrelation [16] and atmospheric delay [17], scholars have proposed timeseries InSAR techniques, such as Persistent Scatterer InSAR (PS-InSAR) [18] and Short Baseline Subset InSAR (SBAS-InSAR) [19]. The PS-InSAR or SBAS-InSAR techniques can recover long-temporal time-series ground deformations with accuracy levels of mm to cm. Recently, scholars have applied time-series InSAR to monitor the landslide at hydropower stations, such as the Three Gorges Dam [10] and the Wudongde Hydropower Station [20].

The Xiluodu Hydropower Station is located upstream of the Jinsha River in China. As the world's third-largest hydropower station, the Xiluodu Hydropower Station has experienced frequent landslides and other geohazards since its impoundment in 2013. In 2013, the Huangping landslide, with an inflow of earthworks amounting to <sup>20</sup> × 104 m3, caused 12 deaths. The Ganhaizi landslide, with an influx of earthworks of 7800 × <sup>104</sup> <sup>m</sup>3, was the largest landslide that occurred at the Xiluodu Hydropower Station [21]. LIANG Guohe [22] and ZHOU Zhifang [23] used the leveling points and valley-shrink surveying lines, respectively, to monitor subsidence and valley-shrink deformation around the shoreside of the Xiluodu reservoir and the Xiluodu arch dam, and their results showed that the Xiluodu arch dam and the shoreside of the reservoir were under valley-shrinkage deformation of 50 mm, accompanied by about 25 mm subsidence. Additionally, ZHOU Zhifang [23] pointed out that valley-shrinkage deformation is highly correlated with the water level. In 2017, LI Lingjing [24] analyzed the Yizicun landslide with ALOS, ASAR, and TerraSAR datasets and found that there was a push-type landslide before the water impoundment and that the landslide changed to a pull-type after the water impoundment.

Currently, the research on the Xiluodu Hydropower Station has focused on valleyshrinkage deformation of the shoreside and the arch dam and their subsidence. In 2019, Li et al. [25] used D-InSAR technology to determine the vertical and horizontal displacement of the landslide near the Xiluodu Hydropower Station and further analyzed its deformation characteristics. In 2022, Zhu et al. [26] obtained the horizontal and vertical displacement of the landslide of the Xiluodu Water Station by SBAS technology, using ALOS2 ascending and descending data, and explained the cause of slope instability by combining geomorphic data and lithologic characteristics.

However, these studies were mainly based on traditional measurements and could not completely reflect the ground deformations occurring near Xiluodu. Although the InSAR technique was used to monitor the Xiluodu reservoir area, it was limited to a specific landslide. Furthermore, the hydrogeological conditions have changed after water impoundment. Whether the change in the hydrogeological conditions will lead to widerange deformation in Yongshan and Leibo counties and change the relationship between the water level and deformation, needs to be analyzed.

#### **2. Geological Setting**

#### *2.1. Geological Setting of Xiluodu*

As the backbone project of "Power Transmission from West to East" in China, Xiluodu is the largest hydropower station on the Jinsha River. The Xiluodu Hydropower Station is located in the U-shaped river valley upstream of the Jinsha River. On its left side is Leibo County of Sichuan Province, and on its right side is Yongshan County of Yunnan Province. The construction of Xiluodu Hydropower Station started in April 2007, and Xiluodu Hydropower Station carried out the first stage of water storage trial operation in May 2013. In June 2014, Xiluodu Hydropower Station officially began operating, and at the end of September, the water level of the reservoir reached 600 m for the first time. The Xiluodu Hydropower Station is mainly designed to generate power, control floods, fix sand, and navigation [27].

As shown in Figure 1, the reservoir area is located on the southwestern margin of the Yangtze platform, which is the transitional zone between the Yunnan–Guizhou Plateau and the Sichuan Basin. The terrain is generally high in the west and low in the east, belonging to the strong erosion of the high-mountain landform [28]. The Xiluodu Hydropower Station is located at the junction of the first and second steps of China. Rapid changes in the terrain may easily cause geological disasters. Moreover, with the completion of the hydropower station, significant fluctuations in water levels will significantly affect the fragile geological environment.

**Figure 1.** Topography profile of Western—Eastern China.

As shown in Figure 2, the Xiluodu area is surrounded by the yellow mud slope anticline in the north, the wall rock syncline in the south, the Majinzi fault (F1) in the west, and the Jiziba fault (F2) in the east, forming a relatively stable zone. However, a thrust—nappe structure fault (named Majiahe, F3) is located 3 km from the dam area, and the number of thrust outliers is influenced by the Majiahe fault. When the water level of the Xiluodu reservoir increases, most of the thrust outliers are submerged, and subsequently, the fault becomes unstable, affecting the safety of the dam area [29].

#### *2.2. Impoundment of Xiluodu*

As shown in Figure 3a, the water level reached its highest level of 600 m in October 2014 for the first time. The water level dropped to 590 m in March 2015. In June 2015, the water level rapidly dropped to 540 m, and then a new round of water storage was started. In July 2015, the water level rose to 555 m. In August 2015, the water level rose to 560 m. In October 2015, it reached 600 m again. Since the water impoundment in October 2014, four complete water impoundment cycles have been completed. Because the annual rainfall is different, the water storage mode is adjusted accordingly.

**Figure 2.** Basic geological structural map of Xiluodu area. F1 is the Majingzi fault, F2 is the Jiziba fault, and F3 is the Majiahe fault.

**Figure 3.** (**a**) Water level of Xiluodu reservoir; (**b**) distribution of aquifers (along the red line in Figure 2) of Xiluodu after water impoundment. A is the river reservoir, B is the area of Yongshan County. S is the Rammell aquifuge layer; P1 m is the limestone permeable bed layer; P2βn is the Rammell aquifuge layer; P2β is the basalt permeable bed layer; P2x is the sandy shale layer; Q is the Quaternary loose accumulation layer; the green arrow indicates the direction of water flow; the red arrow indicates the direction of deformation along the Jinsha River valley. L1 is the water level before constructing Xiluodu Hydropower Station, and L2 is the water level after its construction.

The distribution of the Xiluodu aquifers changed significantly after water impoundment. The geological setting mainly comprises three water-repellent layers (i.e., S, P2βn, and P2x) and three water-permeable layers (i.e., P1 m, P2β, and Q) [23]; see Figure 3b. The highest water level is 600 m, which is lower than the elevation of the sandy shale layer (P2x). Therefore, only the P1 m and P2β permeable bed layers are directly connected with

the reservoir water during the process of impoundment. The potential energy of Xiluodu water increases with the water level. For the P1 m limestone permeable bed layers in a pressurized state, the water permeates the P1 m layer through the exposed limestone area, increasing the P1 m pressure potential energy. For the P2β basalt permeable bed layers in a state of low pressure, the water will directly flow into the P2β layer through the submerged basalt area (see the green arrow in Figure 3b), causing the P2β layer to be saturated gradually.

When the water level rises to 600 m, the rising value of the water level (the difference in L2 and L1) is about 220 m. It leads to an increase in the water's potential energy. Due to the existence of confining beds (Rammell layer P2βn), the groundwater pressure potential energy under P2βn will generate an oblique upward pressure (see the solid red arrow in Figure 3b), which will cause the uplift and horizontal displacement perpendicular to the shoreline (see the red dash arrow in Figure 3b). Moreover, as the P2β layer becomes increasingly saturated, the saturated soil will consolidate under the action of gravity, causing vertical subsidence.

#### **3. Methods**

#### *3.1. Theoretical Basis of Two-Dimensional Deformation Decomposition*

The InSAR measures a one-dimensional deformation along the line of sight (LOS) direction, and such measurements can lead to misinterpretation, as the true deformations may consist of horizontal and vertical displacements [30]. As is known, the LOS deformation is the sum of projection results of horizontal (including north–south and east–west components) and vertical displacement in the LOS direction. Given that there are three or more LOS measurements from different viewing geometries, we can decompose the LOS measurements into north–south, east–west, and vertical displacements by least-squares adjustment [31]. However, for this specific study, we know that displacement along a certain direction is zero or can be ignored based on certain assumptions. For the Xiluodu Reservoir, which is a valley reservoir, the displacement along the direction parallel to the shoreline can be ignored. Thus, the horizontal displacement in a direction perpendicular to the shoreline and the displacement in a vertical direction are sufficient for interpretation. Moreover, due to the small number of viewing angles in the study area, it lacks any redundant observation for us to check for errors if we decomposed the three-dimensional deformations. Thus, we proposed an approach to decompose the InSAR measurements in the LOS direction into two-dimensional displacements, i.e., the horizontal displacement in a direction perpendicular to the shoreline and the vertical displacement.

The specific geometry for decomposing LOS deformation into two-dimensional displacements is shown in Figure 4.

The LOS deformation measurement can be expressed as

$$d\_{l\infty} = d\_{l\text{\textquotedblleft}} \cos \theta + d\_{l'} \cos(\delta - \varphi) \sin \theta \tag{1}$$

where *δ* is the azimuth angle of a predefined direction, i.e., the clockwise angle with respect to the north direction; *θ* is the incident angle of the radar; *ϕ* is the heading angle of the satellite; *h* and *u* are the horizontal deformation along predefined direction and vertical deformation, respectively; and *dlos* is the InSAR measured deformation in the LOS direction.

Given that there are *K* (*K* >= 2) LOS measurements for some specific time periods, there is a set of observation equations written in matrix form:

$$
\begin{bmatrix} d\_{los1} \\ d\_{los2} \\ \vdots \\ d\_{losK} \end{bmatrix} = \begin{bmatrix} \cos(\delta - \varphi\_1)\sin\theta\_1 & \cos\theta\_1 \\ \cos(\delta - \varphi\_2)\sin\theta\_2 & \cos\theta\_2 \\ \vdots & \vdots \\ \cos(\delta - \varphi\_K)\sin\theta\_K & \cos\theta\_K \end{bmatrix} \times \begin{bmatrix} h \\ u \end{bmatrix} \tag{2}
$$

$$Ax = b \tag{3}
$$

where *A* is the design matrix with a size of *K* × 2, *x* is the unknown parameter vector, and b is the observation vector.

When solving Formula (3), because there are random errors in the interferogram, such as noise and some unremoved clean atmospheric delay phase, a random model needs to be added to constrain some of the errors. In this paper, interferogram coherence is used for weighting; then, we can obtain the two-dimensional deformations by using the least-squares adjustment method:

$$\mathbf{x} = \left(A^T P \mathbf{A}\right)^{-1} A^T P \mathbf{b} \tag{4}$$

where *P* is the weight matrix composed of interferogram coherence.

**Figure 4.** Schematic diagram of LOS decomposition. The red arrow represents the local east, north, up (ENU) coordinates; the blue arrow represents the flight direction of the satellite; the green arrow indicates the LOS direction; *ϕ* and *θ* are the azimuth and incident angle of the satellite, respectively; *δ* is the azimuth angle of a predefined direction; *n* and *p* directions are defined local coordinate systems.

#### *3.2. InSAR Two-Dimensional Deformation Method*

Recently, the decomposition of three-dimensional deformation sequences has been developed. Samsonov and D'Oreye proposed the MSBAS method to decompose the timeseries East–West and vertical deformations [32] by ignoring the North–South deformations; Pepe et al. proposed the minimum-acceleration (MinA) method, which is used to recover three-dimensional (3D) time-series deformations [33]; Xuguo Shi et al. recovered threedimensional deformation sequences by using cubic spline interpolation [34]. However, rather than obtaining the 3D deformations or the East–West and vertical deformations by ignoring the North–South deformations, we obtained the vertical and horizontal deformation, called two-dimensional (2D) deformations hereinafter, along a specified direction in reality, which benefits the interpretation.

As shown in Figure 5, given that we have three tracks of SAR data, that is, the Sentinel-1 ascending and descending and ALOS2, covering the study area. The SAR datasets are processed to obtain the time-series LOS deformations. The deformations are geocoded into a geographic coordinate system and further resampled to a common geographic grid. Assuming that the number of images of the three tracks SAR data is *N*1, *N*2, and *N*3, and the acquisition times are *T*<sup>1</sup> = *t* (1) <sup>1</sup> , *t* (1) <sup>2</sup> ,..., *t* (1) *N*1 , *T*<sup>2</sup> = *t* (2) <sup>1</sup> , *t* (2) <sup>2</sup> ,..., *t* (2) *N*<sup>2</sup> , and *T*<sup>3</sup> = *t* (3) <sup>1</sup> , *t* (3) <sup>2</sup> ,..., *t* (3) *N*<sup>3</sup> , respectively.

**Figure 5.** Time series of Sentinel-1 ascending and descending data, and ALOS2 ascending data.

Assuming that the acquisition times are unique for different pieces of SAR data, the total number of images is *N* = *N*<sup>1</sup> + *N*<sup>2</sup> + *N*3. The times of all of the data acquisitions are *T* = *T*<sup>1</sup> ∪ *T*<sup>2</sup> ∪ *T*3. The corresponding LOS time-series deformations sequence are *d*<sup>1</sup> = *d* (1) <sup>1</sup> , *d* (1) <sup>2</sup> ,..., *d* (1) *N*1 , *d*<sup>2</sup> = *d* (2) <sup>1</sup> , *d* (2) <sup>2</sup> ,..., *d* (2) *N*<sup>2</sup> , and *d*<sup>3</sup> = *d* (3) <sup>1</sup> , *d* (3) <sup>2</sup> ,..., *d* (3) *N*<sup>3</sup> ; then, the merged LOS time-series deformations are *d* = *d*<sup>1</sup> ∪ *d*<sup>2</sup> ∪ *d*3.

According to Equation (2), the LOS time-series deformations d can be written as a set of equations with respect to the time-series 2-D deformations, and expressed in matrix form:

$$G\mathfrak{x} = d\tag{5}$$

where *d* = *d*<sup>1</sup> *d*<sup>2</sup> ··· *dN <sup>T</sup>* is the observation vector, and the unknown parameters vector *x* = *h*<sup>1</sup> *u*<sup>1</sup> *h*<sup>2</sup> *u*<sup>2</sup> ··· *hN uN <sup>T</sup>* with a size of 2*<sup>N</sup>* <sup>×</sup> 1. The design matrix has a size of *N* × 2*N*, and its form is as follows:

$$G = \begin{bmatrix} A\_1 \\ & A\_2 \\ & & A\_3 \\ & & & \ddots \\ & & & & A\_N \end{bmatrix}$$

where *A*<sup>1</sup> = *cos*(*<sup>δ</sup>* <sup>−</sup> *<sup>ϕ</sup>*(1))*sinθ*(1) *cosθ*(1) , *A*<sup>2</sup> = *cos*(*<sup>δ</sup>* <sup>−</sup> *<sup>ϕ</sup>*(2))*sinθ*(2) *cosθ*(2) , and *A*<sup>3</sup> = *cos*(*<sup>δ</sup>* <sup>−</sup> *<sup>ϕ</sup>*(3))*sinθ*(3) *cosθ*(3) . The superscript in the angles of *ϕ* and *θ* is dependent on which subset of *di* it is from.

Obviously, Equation (5) is rank deficient as the number of equations is smaller than the number of parameters to be estimated, i.e., *N* < 2*N*. To solve the problem of rank deficiency, we add smoothing constraints to the parameters by assuming that the deformation rates of adjacent displacements are equal, i.e.,

$$\frac{h\_{i+2} - h\_{i+1}}{t\_{i+2} - t\_{i+1}} = \frac{h\_{i+1} - h\_i}{t\_{i+1} - t\_i} \tag{6}$$

$$\frac{u\_{i+2} - u\_{i+1}}{t\_{i+2} - t\_{i+1}} = \frac{u\_{i+1} - u\_i}{t\_{i+1} - t\_i} \tag{7}$$

where *i* ∈ [1, *N* − 2]. Then, we obtain a set of equations with respect to unknown parameters vector, and it is expressed by

$$
\mathbb{R}\mathfrak{x} = 0\tag{8}
$$

where *R* is the roughness matrix with a size of (*N* − 2) × *N*, and


Combining Equations (5) and (8), we solve the regularized least-squares problem [35]:

$$\min \left\| \begin{bmatrix} G \\ \alpha \mathcal{R} \end{bmatrix} x - \begin{bmatrix} d \\ 0 \end{bmatrix} \right\|\_{2}^{2} \tag{9}$$

The Tikhonov regularization solution is determined by

$$\mathbf{x} = \left(\mathbf{G}^T \mathbf{G} + \boldsymbol{\alpha}^2 \mathbf{R}^T \mathbf{R}\right)^{-1} \mathbf{G}^T \mathbf{d} \tag{10}$$

where *α* is the smoothing parameter and is unknown beforehand. We determine it by using the L-curve method [35].

Figure 6 presents a flowchart of the proposed method.

**Figure 6.** Flowchart of the proposed method.

#### **4. Data Processing and Results**

#### *4.1. SAR Data Used*

In order to study surface deformation near the Xiluodu reservoir area, we used the ascending and descending data of Sentinel-1 and the ascending data of ALOS2: from October 2014 to March 2019, there was a total of 100 ascending Sentinel-1 data, and the incident angle of the image center was 39.29◦. From October 2014 to April 2019, there were 97 descending Sentinel-1 data, and the image center incident angle was 33.9◦. From September 2014 to December 2018, there was a total of 12 pieces of ascending ALOS2 data, and the image center incident angle was 31.4◦; the coverage area is shown in Figure 7. See Table 1 for specific information. At the same time, we chose the Shuttle Radar Topography Mission (SRTM) with a 30 m x 30 m resolution to remove the terrain phase.

**Figure 7.** (**a**) Study area; (**b**) SAR data coverage area. The white dashed box represents the Sentinel-1 ascending data, and the black dashed box represents the Sentinel-1 descending data. The blue dashed box represents the ALOS2 ascending data.


**Table 1.** Three kinds of data for SBAS interference pair information.

#### *4.2. Data Processing*

We first performed SBAS processing on each piece of data to obtain the deformation sequence of Sentinel-1 data and ALOS2 data. In the data-processing phase, first, in order to ensure the consistency of the three sets of data results, we performed 8 × 2 multi-look on the Sentinel-1 data and 3 × 7 multi-look on the ALOS2 data, thus obtaining image results of about 15 m resolution. Second, because the Sentinel-1 data and ALOS2 data orbital control were accurate, the spatial baseline was short, so the influence of spatial decoherence was not considered. At the same time, in order to overcome the time decoherence, we adopted the two-connected interference pairing method (as shown in Figure 8: 1 and 2 constitute the interference pair, 1 and 3 constitute the interference pair, 2 and 3 constitute the interference pair, and so on), so that any scene image was in a triangular closed loop, which overcame the time decoherence and could test for errors. Therefore, we obtained 21 interference pairs from the ALOS2 ascending data, 197 interference pairs from the Sentinel-1 ascending data obtained, and 191 interference pairs from the Sentinel-1 descending data. Third, we used the SRTM data of 30 × 30 m resolution to obtain the DEM in the SAR coordinate system by geocoding, thus obtaining the simulated terrain phase, and then using the generated interferogram subtracts in the simulated interference phase to remove the terrain phase. In order to eliminate the phase residual and improve the signal-to-noise ratio of the interferogram, we used a 'Goldstein' filter on the differential interference phase with a filter factor of 0.5 and a filter window size of 64. As the low-coherence region could be affected by phase noise, we averaged the coherence coefficient map when unwrapping and then ensured that the region with coherence below 0.52 did not participate in the unwrapping, thus ensuring the reliability of the wrap. In the phase unwrapping stage, we selected the unwrapping reference point in the stable region away from the city. After unwrapping with the minimum cost flow, we refined the baseline of the distracted differential interferogram, thereby effectively removing the track phase. Through the above process, we obtained the exact differential interference phase after unwrapping and then converted it into a settlement sequence result. Finally, we inversely geocoded the resulting sedimentation results and resampled them to the same geographic coordinate range to obtain the deformation sequence results for the three images at 15 m resolution.

**Figure 8.** Two-connection network mode. The black triangles represent the SAR images in chronological order from left to right, and the red lines represent different interference pairs.

On this basis, using the algorithm described in Section 3.2, the deformation sequence was decomposed to obtain the deformation sequence results in the vertical shore direction and the vertical direction, and the deformation rate in the horizontal direction and the vertical direction was calculated by using the two-dimensional deformation sequence. The results are provided in Section 4.3.

#### *4.3. Results*

#### 4.3.1. LOS Deformation

Figure 9 shows the average rate map of the Sentinel-1 data and ALOS2 data obtained using the SBAS technique. The reference point was selected in the same stable area farther from the deformation zone. For the Yongshan County area, the ALOS2 ascending data LOS has a maximum average speed of −90 mm/a; the ascending Sentinel-1 data LOS has a maximum average speed of −110 mm/a; and the descending Sentinel-1 data LOS has a maximum average rate of 60 mm/a (positive values represent proximity to the satellite and negative values represent the distance from the satellite). For the ALOS2 ascending data and the Sentinel-1 ascending data, the incident angles are 31.4◦ and 39.29◦, respectively. The reasons for the different deformation rates may be as follows: (1) the data incident angles are different, resulting in different observational geometry; (2) Xiluodu horizontal shifts may occur in the reservoir area, resulting in different projections onto the LOS backward direction. For the Sentinel-1 ascending data and the descending Sentinel-1 data, the deformation rate results are opposite, indicating that the horizontal displacement

must occur in the Yongshan County area. Therefore, we needed to achieve this in order to decompose and obtain horizontal displacement and vertical displacement.

**Figure 9.** (**a**) Mean velocity map of ALOS2 ascending; (**b**) mean velocity map of Sentinel-1 ascending; (**c**) mean velocity map of Sentinel-1 descending.

#### 4.3.2. Two-Dimensional Deformation Rate

We used the Section 3.2 formula to decompose the deformation sequence to obtain the deformation sequence results in the horizontal direction and the vertical direction and then obtained the mean velocity map perpendicular to the shore direction and the vertical direction, as shown in Figure 10. As can be seen from Figure 10, the maximum deformation rate is 250 mm/a perpendicular to the outward direction of the slope, and the maximum deformation rate is 60 mm/a in the vertical downward direction. It can be seen that the deformation of the Xiluodu reservoir area is mainly in the horizontal direction, accompanied by vertical sedimentation, which is consistent with the previous results of Zhou et al. [23]. Furthermore, we sought to determine the relationship between water storage and horizontal and vertical directions by analyzing the deformation sequence.

**Figure 10.** (**a**) Mean rate perpendicular to the shore direction; (**b**) vertical mean rate.

In order to more intuitively display the surface deformation of the Xiluodu reservoir area, because the number of images was too large, we drew some horizontal and vertical deformation sequence images, as shown in Figure 10. It can be seen from Figure 11 that, since October 2014, the horizontal deformation of the Yongshan County area has been accelerating, and the maximum deformation in April 2019 was about 1000 mm; at the same time, the vertical direction deformation has also slowly increased, and by April 2019, the

maximum deformation was about 200 mm. Since the reservoir was filled with water in 2014, such a huge deformation occurred in the vicinity of the Xiluodu reservoir area. According to previous studies, the reservoir impoundment may lead to instability of the reservoir slope [36]. Therefore, we suspect that due to the huge water level after the impoundment, the water level difference caused by the change threatens the Xiluodu surface, which could lead to the deformation of the Earth's surface.

**Figure 11.** (**a**) Horizontal time-series deformation; (**b**) vertical time-series deformation.

#### 4.3.3. Two-Dimensional Time Series

It can be seen from Figure 10 that the closer the bank is, the smaller the deformation rate, and the farther away it is from the shore, the greater the deformation rate. At the same time, the study area is obviously two distinct regions. In order to prevent the formation of the trailing edge of the study area from creating a landslide and to prevent the formation of a ground crack in the middle of the study area, we selected the section line AB, CD perpendicular to the shore direction and the section line EF parallel to the shore for analysis. Figure 12a–d are the deformation rate profile lines of the vertical shore, and Figure 12e,f are the deformation rate profile lines of the parallel shore. (1) The AB section line reaches a peak at 3.7 km from the shore, while the horizontal deformation rate is greater than 100 mm and the range is about 0.75 km, and the vertical deformation rate is greater than −40 mm. The range is about 0.5 km. (2) The CD profile line reaches a peak at 3.4 km from the shore, while the horizontal deformation rate is greater than 100 mm, which is about 0.65 km; the vertical deformation rate is greater than −20 mm, which is about 0.6 km. (3) The EF section achieved peaks at 1.8 km and 4.2 km from point E. It can be seen from the results that these two areas are more dangerous, and we will continue to analyze them in combination with optical images.

**Figure 12.** (**a**) is the horizontal displacement rate of profile line A–B; (**b**) is the vertical displacement rate of profile line A–B; (**c**) is the horizontal displacement rate of profile line C–D; (**d**) is the vertical displacement rate of profile line C–D; (**e**) is the horizontal displacement rate of profile line E–F; (**f**) is the vertical displacement rate of profile line E–F. The green line represents the suspected fault.

#### **5. Analysis and Discussion**

#### *5.1. Precision Analysis*

The accuracy of this experiment was evaluated by calculating the RMS of the residual of the deformation sequence. The results are shown in Figure 13. Furthermore, we statistically analyzed the RMS of the deformation sequence (see Table 2 for details). We found that the mean root mean square error of the region is 4.7 mm, RMS < 10 mm is about 95.36% of the total, and RMS < 15 mm is about 98.85% of the total. The root error average is 4.1 mm, RMS < 10 mm accounts for 95.40% of the total, and RMS < 15 mm accounts for 99.31% of the total. The mean square root error of the non-deformation significant region is 4.8 mm, RMS < 10 mm accounts for 95.35% of the total, and RMS < 15 mm is approximately 98.79% of the total. It can be seen from Figure 14 that the deformation region mainly occurs in the county area, and the terrain is relatively flat, so the accuracy of the InSAR monitoring result

is high, while the other regions have steep terrain and relatively low precision. Therefore, we think that the deformation accuracy could reach 10 mm within the 95% confidence interval, so the accuracy of this experiment is reliable.

**Figure 13.** Root mean square error chart. (**a**) Total mean square error; (**b**) deformation significant regional mean square error; (**c**) deformation non-significant area mean square error.

**Table 2.** RMS data statistics.


**Figure 14.** Overall root mean square error map.

#### *5.2. Relationships between Displacement and Water Storage*

According to the literature [23], the elevation of the Jinsha River riverbed near the Xiluodu reservoir area was about 400 m before water storage. When the Xiluodu Hydropower Station began to store water, the lowest water level in the reservoir area rose to 540 m, and the highest water level rose to 600 m. It leads to a continuous increase in surface deformation near the reservoir area. At the same time, the fluctuation of the Xiluodu water level of about 60 m per year will also cause the surface deformation near the reservoir area

to fluctuate accordingly. To this end, we selected two feature points, P1 and P2, in the two significant areas of deformation, to analyze the relationship between water storage and deformation and provide a basis for future reservoir water storage.

#### 5.2.1. Relationship between Horizontal Displacement and Water Storage

As shown in Figure 15, the water of Xiluodu has been stored for the first time since October 2014, and the water level dropped to 540 m in June of the following year. At present, it has experienced four complete water storage–drainage cycles, of which the cumulative shape variable of the direction of the P1 level reached 110 mm, and the cumulative shape variable of the horizontal direction of the P2 level reached 800 mm. At the same time, the horizontal deformation rate obviously accelerated in each round of water storage, and when the water level of Xiluodu rises to 600 m, the horizontal deformation rate begins to decrease. Over time, in the third round and the fourth round of water storage, although the absolute deformation variable increased in the horizontal direction of the Xiluodu reservoir area, the relative deformation rate began to slow down, indicating that the horizontal deformation of the Xiluodu reservoir area tends to be stable. According to Zhou et al.'s literature analysis [23], as the water level continues to rise, the pore water pressure inside the permeable layer of Yangxin limestone increases continuously, which leads to a reduction in effective stress, which leads to the elastic expansion deformation of Yangxin limestone. At the same time, since the upper layer of the Yangxin limestone is covered with mud shale as a water-blocking layer, the shale also undergoes elastic deformation as the pressure increases. Since the Yongsheng syncline basin generally tends toward the Jinsha River, the horizontal direction of the Xiluodu reservoir area will be deformed outwardly along the slope under the superposition of two elastic deformations. However, as the water level of Xiluodu continues to rise, the water level of the reservoir enters from the permeable layer of basalt, resulting in a stable horizontal deformation.

**Figure 15.** Relationship between horizontal displacement and impoundment. (**a**) is the relationship between the deformation sequence and water level in the horizontal direction of point P1 in Figure 10; (**b**) is the deformation sequence in (**a**) that removes the linear trend; (**c**) is the relationship between the deformation sequence and water level in the horizontal direction of point P2 in Figure 10; (**d**) is the deformation sequence in (**c**) that removes the linear trend.

#### 5.2.2. Relationship between Vertical Displacement and Water Storage

As shown in Figure 16, from October 2014 to April 2019, the cumulative shape variables at the P1 point and the vertical point of the P2 point were −175 mm and −110 mm, respectively. At the same time, when the water level of the reservoir rapidly decreases, the deformation in the vertical direction accelerates. When the water level of the reservoir rises, the deformation in the vertical direction begins to rebound and finally, the deformation in the vertical direction changes with the fluctuation of the water level. Therefore, we speculate that when the water in the Jinsha River continues to infiltrate, the floating force is generated, causing the Xiluodu reservoir area to move upward in a vertical direction. However, the elevation of the water level leads to an increase in the gravitational potential energy, and the soil is consolidated under the action of gravity. The gravity is much greater than the buoyancy. Therefore, the Xiluodu reservoir area has a vertical downward deformation. Finally, as the water-storage cycle continues to increase, the natural consolidation state of the soil becomes saturated so that the cyclical changes caused by the water storage are not offset; thus, the vertical deformation changes with the water level fluctuation.

**Figure 16.** Relationship between vertical displacement and water storage. (**a**) is the relationship between the deformation sequence and water level in the vertical direction of point P1 in Figure 10; (**b**) is the deformation sequence in (**a**) that removes the linear trend; (**c**) is the relationship between the deformation sequence and water level in the vertical direction of point P2 in Figure 10; (**d**) is the deformation sequence in (**c**) that removes the linear trend.

#### *5.3. Potential Geological Hazards*

According to the National Geological Hazard Bulletin, there were 2966 geological disasters in 2018, resulting in 105 deaths and seven missing people. The direct economic losses amounted to RMB 1.47 billion. Among them, the southwestern region suffered the most serious disasters, accounting for 35.2% of national geological disasters. Therefore, monitoring geological disasters is vital. As a wide-ranging monitoring method, InSAR can undoubtedly provide necessary observation of geological disasters. However, due to the limitation of satellite side-view imaging, the SAR satellite of a single platform may have overlapping, shadow, and top–bottom inversion in complex mountainous areas, which

may easily cause blind spots of identification, resulting in misjudgment [37]. Therefore, through a variety of data (ascending and descending data), the true strength of InSAR technology can be realized in complex mountainous areas. In this paper, we used the data of multiple orbits to decompose the study area horizontally and vertically, thus restoring the deformation of the entire study area. According to the analysis in Section 5.1, combined with the optical image (as shown in Figure 17a), we found that the deformation of the trailing edge (fault-1) of the significant region of deformation was large, and then through the analysis in Section 5.2, we speculate that the region will continue to move along the coast, which may lead to cracks at the trailing edge, and the formation of geological disasters such as ground fissures and collapses, which require following up. At the same time, the left side of the area has a distinct block in the deformed area on the right. As the deformation in the horizontal direction continues to accumulate, the unequal horizontal forces in the two areas may cause relative shear, causing the middle area to form a fault (Fault-2), which may lead to geological disasters such as landslides and ground fissures. As can be seen from the foregoing, the geological conditions in the area are relatively fragile. Therefore, in order to avoid large casualties, the area should be continuously monitored.

**Figure 17.** (**a**) Optical image; (**b**) potential fault.

#### **6. Conclusions**

As China's demand for power resources continues to increase, hydropower has become increasingly important as a clean and recyclable energy source. However, the artificial construction of a reservoir undoubtedly changes the hydrogeological conditions of the hydropower station's attached structures, resulting in deformation of the surface near the hydropower station, which can lead to landslides, mudslides, earthquakes, etc. In this paper, using time-series InSAR technology, through a variety of data (ALOS2 ascending, Sentinel-1 ascending and descending), large-scale surface deformation monitoring was carried out near the Xiluodu reservoir area after water storage and power generation. The conclusions are as follows:

(1) For a variety of data, the space is not synchronized, and the deformation sequence decomposition cannot be directly performed. A new solution was proposed, and the horizontal and vertical deformation values of the entire time series were obtained. The results show that the maximum average velocity is perpendicular to the bank edge, the maximum average velocity is 250 mm/a, and the maximum vertical velocity is 60 mm/a, which indicates that the surface deformation near the Xiluodu reservoir area is mainly horizontal and the vertical direction is supplemented. At the same time, we used the residuals to analyze the accuracy of the deformation solution. The results show that the deformation accuracy can reach 10 mm within the 95% confidence interval, indicating that the accuracy of this experiment is reliable.


This experiment failed to obtain GPS data, so it is impossible to verify the decomposition results of InSAR data with external data. At the same time, the acquisition date of this experimental data was after the Xiluodu Hydropower Station started operating. In subsequent analysis, we can try to obtain the InSAR data before the Xiluodu Reservoir was stored and thus gain a clearer understanding of the landmark deformation near the Xiluodu Reservoir.

**Author Contributions:** Conceptualization, Q.C., B.X. and Z.L.; data curation, H.Z.; formal analysis, W.M. and Z.L.; funding acquisition, B.X.; investigation, W.M., Z.L. and B.X.; methodology, Z.L. and Q.C.; supervision, B.X.; writing—original draft, Q.C. and Z.L.; writing—review and editing, all authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partly supported by the National Science Fund for Distinguished Young Scholars (No. 41925016), and the National Natural Science Foundation of China (No. 41804008).

**Data Availability Statement:** The Sentinel-1 data have been made freely available by the European Space Agency and distributed and archived by the Alaska Satellite Facility (https://www.asf.alaska. edu/sentinel/, 25 October 2019). The ALOS2 SAR data were provided by JAXA under the second Research Announcement on the Earth Observations (EO-RA2) for its Earth observation satellite projects (Project No. ER2 A2 N167).

**Acknowledgments:** Several figures were plotted by General Mapping Tools (GMT v6.0).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Yafei Wang 1,2, Yao He 1,2, Jiuyi Li 1,3 and Yazhen Jiang 1,4,\***


**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<sup>2</sup> 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.

**Keywords:** GEE; long-term time series images; land use classification; spatial modeling; land degradation; Huangshui watershed

#### **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 highlevel 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

**Citation:** Wang, Y.; He, Y.; Li, J.; Jiang, Y. Evolution Simulation and Risk Analysis of Land Use Functions and Structures in Ecologically Fragile Watersheds. *Remote Sens.* **2022**, *14*, 5521. https://doi.org/10.3390/ rs14215521

Academic Editors: Qiusheng Wu, Xinyi Shen, Jun Li and Chengye Zhang

Received: 29 September 2022 Accepted: 31 October 2022 Published: 2 November 2022

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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 longterm 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 fullelement 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/km2, 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 km2 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.

**Figure 1.** Overview map of the study area (WGS84, 1:1,000,000).

#### *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.

**Figure 2.** The main technical flows adopted in the methodology.

#### *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 longterm 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 *wm* and the seed window *ws* 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 *wm* immediately outside the seed window *ws*; (3) determining the dominant classification *fm* in the sliding window *wm* by investing the frequency  of a function classification: if it was >0.7, the function classification was regarded as the dominant classification *fm* of the sliding window *wm*, 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 *ws* was consistent with the dominant classification *fm* in the adjacent sliding detection window *wm*: if they were consistent, we set all the classifications in *wm* to *fs*, and moved the position of the seed window to the grid position where the last dominant classification *fs* 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 *Pm*,*<sup>q</sup> <sup>t</sup>* as a function parameter year by year. The *Pm*,*<sup>q</sup> <sup>t</sup>* was expressed as follows:

$$P\_t^{m\mathcal{A}} = \frac{s\_t^{m\mathcal{A}}}{s}$$

where *s m*,*q <sup>t</sup>* 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−<sup>4</sup> 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 *Sk* of land-use. The functional suitability probability *Sk* of land-use was expressed as follows:

$$S\_k = \sum\_{\mathbf{x}} w\_{\mathbf{x},k} \star I\_{\mathbf{x}}$$

where *x* is the impact factor, *wx*,*<sup>k</sup>* is the impact weight of the impact factor *x* on the functional classification *k*, and *Ix* is the quantitative parameter of the impact factor *x*.

Using *Sk* as the driving factor, and the land scale simulation result of functional classification as the target parameter, the generation probability *TP<sup>t</sup> <sup>p</sup>*,*<sup>k</sup>* 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].

$$TP\_{p,k}^t = \mathcal{S}\_{p,k} \* \mathcal{Q}\_{p,k}^t \* Intritai\_k^t$$

where *k* is the land type; *p* is the spatial location; *t* represents the *t*-th iteration; *TP<sup>t</sup> <sup>p</sup>*,*<sup>k</sup>* represents the probability that the functional classification; *Sp*,*<sup>k</sup>* is the functional suitability probability of land-use; *Ω<sup>t</sup> <sup>p</sup>*,*<sup>k</sup>* represents the neighborhood influence; and *Intertia<sup>t</sup> <sup>k</sup>* 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:

$$\begin{cases} Z = f(X, Y) \\ X = \mathbb{C}\_x(\mathfrak{x}\_1, \mathfrak{x}\_2, \mathfrak{x}\_{3'}, \dots, \mathfrak{x}\_n) \\ Y = \mathbb{C}\_y(y\_1, y\_2, \dots, y\_n) \end{cases}$$

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~*xn* represent different types of ecological service functions, and *y*1~*y*<sup>n</sup> 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 km2. 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 km2 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.

**Figure 3.** Results of year-by-year land use function classification from 1990–2019.

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 km2, accounting for 73.70% of the total area of the watershed; for construction land dominated by human activities as well, it mainly (725.88 km2, 62.72%) located in Xining and Haidong in midstream area. The bare land was mainly distributed upstream, with an area of 348.47 km2, 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 km2 to 1091.83 km2, with an expansion of 2.67 times; while the area of farmland decreased in a wide range, from 6084.17 km2 to 4900.33 km2, 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 km2, 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 km2, 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<sup>2</sup> and 160.21 km2 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. 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 km2, which increased by 330.31 km2 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<sup>2</sup> to 752.20 km2, with an increase of 23.37%. Shrubs and water showed a slight expansion trend, expanding by 277.76 km<sup>2</sup> and 36.13 km2, 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 km2, 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<sup>2</sup> in 2019 to 455.31 km2, with a decline rate of 19.91%. The grassland did not change significantly, and the area change was only 87.07 km2, accounting for 0.82% of the total area.

**Table 1.** Simulation scale result table (km2) + specific gravity.


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<sup>2</sup> and 2938.26 km<sup>2</sup> 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<sup>2</sup> and 121.33 km2, 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 km2, 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 km2, 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 km2, 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 km2 and an expansion rate of 54.54%.

**Figure 4.** Simulation results of land use function and structure in Huangshui watershed in 2035.(WGS84, 1:1,000,000).

#### *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 km2, with a decrease of 4.66%; and the broken farmland with an area of <1000 m2 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.

**Figure 5.** Results of risk analysis of farmland loss: (**a**) Huangshui Basin; (**b**) Duoba Town in Huangzhong District; (**c**) Hetaozhuang town in Minhe County. (WGS84, 1:975,000).

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 km2, 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 km2 (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 km2, 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 km2 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.

**Figure 6.** Results of ecological protection importance evaluation and risk analysis: (**a**) Huangshui Basin; (**b**) Ledu District; (**c**) Lierbao town in Minhe County. (WGS84, 1:975,000).

#### **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–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<sup>2</sup> in 2035, with an increase of 15.73 km2 compared to that in 2019; among them, 14.01 km<sup>2</sup> 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 km2, 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.

**Figure 7.** Distribution of validation sites for risk analysis: (**a**) Duoba Town in Huangzhong District; (**b**) urban area of Ledu District; (**c**) Hetaozhuang town in Minhe County; (**d**) Donggou Town in Huzhu County. (WGS84, 1:975,000).

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 km2, 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 km2, and the simulated degradation area was 9.07 km2, 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 lowcoverage 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.

**Author Contributions:** Conceptualization, Y.W. and Y.J.; formal analysis, Y.W. and Y.H.; funding acquisition, Y.W.; methodology, Y.W. and Y.H.; writing—original draft, Y.W., Y.H., J.L. and Y.J.; writing—review and editing, Y.W., Y.H. and Y.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Second Tibetan Plateau Scientific Expedition and Research Pro-gram (2019QZKK0406), and the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA20020301).

**Data Availability Statement:** Codes available for land cover/land use classification: https://code. earthengine.google.com/674cec07c398a7965fda4e0b1e8df58c (accessed on 7 July 2021). Other parts of the code or process is available from corresponding authors upon reasonable request.

**Acknowledgments:** We appreciate the critical and constructive comments and suggestions from the reviewers that helped improve the quality of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

