Response of Guobu Slope Displacement to Rainfall and Reservoir Water Level with Time-Series InSAR and Wavelet Analysis

Abstract

Reservoir bank landslides are a frequent phenomenon, and the stability of these landslides is affected by two essential factors: rainfall and reservoir level changes. Studying the response patterns of reservoir bank landslide movements to these variables is crucial in preventing their occurrence and mitigating their effects. To this end, this study employed 103 European Space Agency (ESA) Copernicus Sentinel-1 images and the SBAS-InSAR (small baseline subset interferometric synthetic aperture radar) technique to obtain a time series of the Guobu slope deformation from September 2015 to December 2019. The Guobu slope showed significant toppling damage. The satellite line of sight (LOS) detected a maximum subsidence rate of −447 mm/y (the negative sign indicates movement away from the satellite, i.e., subsidence) in the upper section of the slope. Subsequently, three wavelet tools were used to quantitatively analyze the effect of rainfall and reservoir water level on the deformation of the Guobu slope. The results demonstrate a positive correlation between rainfall and the deformation of the Guobu slope. Moreover, the deformation lags behind the rainfall by approximately 70 days. In contrast, the reservoir water level and the deformation of the Guobu slope exhibit an inverse relationship. The deformation of the leading edge of the slope body lags behind the reservoir level by approximately 19 days, while the middle and upper sections of the slope body, which have the most significant rate of variability, lag by about 80 days. Among these factors, rainfall plays a dominant role in the deformation of the Guobu slope, while reservoir levels play a synergistic role. The findings of this study highlight the importance of monitoring and understanding the impact of changes in rainfall and reservoir water levels on the stability of reservoir bank landslides. This understanding is crucial in preventing the occurrence of such landslides and minimizing their impact. The use of remote sensing techniques, together with wavelet analysis, enables the accurate and timely monitoring of the deformation of the Guobu slope, providing valuable insights for disaster warnings and disaster prevention and reduction efforts.

1. Introduction

Reservoir operations (storage and release of water) not only activate already existing paleo landslides, but also induce new landslides, making the reservoir area a geological hazard-prone and high-risk area [1,2,3]. For example, more than 500 landslides occurred between 1941 and 1953 in Roosevelt Lake, resulting from the Grand Coulee Dam construction (USA), causing tremendous economic losses [4]. In October 1963, a landslide occurred on the left bank of the Vaiont reservoir project in Italy; about 2.7 × 10⁸ m³ of the mountain slid into the reservoir, generating a surge of almost 100 m that crossed the dam and destroyed the village and town downstream, causing about 2600 casualties [5,6]. The same problem can be observed in China’s Three Gorges Reservoir Area (TGRA). Since the Three Gorges Reservoir began filling in 2003, the reservoir water level has been raised from 69 m to 175 m ASL (above sea level), and the water level has fluctuated annually between about 145 and 175 m ASL from 2008 to the present, during which a large number of landslides have been induced, such as the Qianjiangping landslide [7], the Shuping landslide [8] and the Muyubao landslide [9]. There have been 4429 geological hazards in the TGRA from Yichang to Jiangjin along the Yangtze River. Among them, 4256 were landslides and rock avalanches, with a total volume of about 4.24 billion m³ [10]. Meanwhile, many scholars have found good agreement between the slow movement of landslides and seasonal rainfall [11,12,13,14]. Therefore, detecting and monitoring the moving bank slopes, evaluating the magnitude of bank slope movement, and analyzing the internal relationship between landslide movement and hydrological factors are of great significance to ensuring the safe operation of hydropower projects.

Landslide deformation is the most intuitive and effective means of characterization of the evolution of landslides, and reveals the inherent risk of landslides [12]. Therefore, monitoring surface deformation can provide critical information for detecting potential geological hazards such as landslides. Interferometry Synthetic Aperture Radar (InSAR) is an emerging technology for regional surface monitoring. Compared with traditional single-point or small-scale landslide deformation monitoring technology, it has the advantages of high monitoring accuracy, wide monitoring range, high spatial resolution, all-day and all-weather operability, etc., and is widely used in landslide research [15,16,17,18,19,20,21]. For example, Grief et al. [15] processed ENVISAT data covering the town of Lubietova in central Slovakia using the Persistent Scatterer (PS) InSAR technique, analyzed the catastrophic landslides that occurred in the area in 1997, determined the precise boundary of activity, and updated the landslide inventory of the area. Liu et al. [16] identified two landslides in Badong County using time-series InSAR technology and found an apparent correlation between the seasonal movement of landslides and water level variation. Yao et al. [17] used small baseline subset (SBAS) InSAR and PS-InSAR techniques to identify 14 landslides in the Jinsha River basin of Gongjue County, and introduced rainfall data to further analyze the spatial–temporal deformation characteristics of two large landslides. Mishra et al. [19] used PS-InSAR technology to delineate and monitor landslide information along the Baglihar Dam reservoir area, located in the Jammu and Kashmir, India. The SAR data consisted of 46 ascending and 70 descending track data acquired from Sentinel-1, and the qualitative analysis of the combined monitoring results was carried out in conjunction with the local precipitation data. All of the studies mentioned above have proven the effectiveness of using InSAR technology in landslide identification and landslide monitoring, and introduced external data (such as earthquake, reservoir water level, rainfall data, etc.) combined with InSAR time series to further analyze the mechanisms of landslide deformation. However, most of them have only qualitatively demonstrated the influence of external factors on landslides, and InSAR monitoring data were not sufficiently utilized. Since wavelet analysis can extend one-dimensional time series to the time and frequency domains, allowing studies to be conducted at multiple scales, this has been widely applied in geoscience [22,23], meteorology [24], praxeology [25], and economics [26] in recent years.

In the first step of this paper, we processed 103 Sentinel-1 images obtained from ESA through the SBAS-InSAR technique to extract information on surface displacements, from September 2015 to December 2019, on the right bank slope, known as the Guobu slope, upstream of the Laxiwa hydropower station. Subsequently, we analyzed the temporal and spatial features of Guobu slope deformation. In the second step, we integrated the InSAR time series results with wavelet analysis to quantitatively investigate the uncertain relationships between rainfall, reservoir levels, and deformed time series data, including the periodicity and seasonality of each time series, along with their correlations and lags. This study elucidates the association between landslide displacement and external factors from a fresh perspective, and its findings can play a vital role in disaster prevention and control.

2. Study Area and Materials

2.1. Study Area

The upper reaches of the Yellow River lie in the transition zone between China’s first- and second-step topographies (Figure 1a). The abundance of water, and the great difference in the topography of the country, mean the water contains high levels of energy. The valley is mainly composed of hard rock masses, including granite and shallow metamorphic rock, with excellent geological engineering properties for dam construction [27]. As a result, more than 20 hydropower stations have been built in the region to meet China’s dual carbon targets and the growing demand for electricity. The Laxiwa hydropower station is the second cascade power station between Longyang Gorge and Qingtong Gorge, located between Guide County and Guinan County in Qinghai, about 32 km from the outlet of the Longyang Gorge Reservoir (Figure 1). The height of the arch dam of the Laxiwa hydropower station is 250 m, the thickness of the dam top is 10 m, the maximum reservoir capacity is about 1 billion m³, the total installed capacity is 4200 MW, and the maximum generation flow is 1900 m³/s. It has the largest installed capacity, the highest generation efficiency, and the largest generating capacity among the numerous hydropower stations in the upper reaches of the Yellow River. The dataset was provided by the National Cryosphere Desert Data Center (http://www.ncdc.ac.cn, access date: 10 August 2021).

The Laxiwa hydropower station was constructed in 2001 and became operational in May 2009 [28]. The first impoundment was completed in March 2009, and the water level increased from 2250 to 2340 m ASL. At the end of May of the same year, the staff found that the Guobu slope had begun to undergo deformation again. The monitoring work began in August of the same year, and 17 GPS monitoring points were arranged on the top platform and surface of the slope; Figure 1e shows their distribution [29]. After an evaluation study was undertaken by relevant experts, the reservoir water level was allowed to rise gradually. Subsequently, the water level was successively raised to 2400 m, 2430 m, and 2448 m ASL in February 2010, February 2011, and September 2012, respectively. The water level remained stable at 2448 m ASL from September 2012 to May 2015, and then was allowed to rise to the designated level of 2452 m ASL at the end of October 2015, and it remained relatively stable after that [30]. Since the first impoundment, the river channels have widened from ~70 m to ~500 m (Figure 1c,d), causing three slopes, including the Guobu slope, to start to move [31]. From 2015 to 2018, two slopes about 4 km and 8 km upstream of the Laxiwa dam exhibited a moderate movement rate of 40–60 mm/y. However, the Guobu slope, showing larger and faster movement, and being situated only 1 km southwest of the Laxiwa hydropower station, has emerged as an imminent threat to the dam, its powerhouse, and the surrounding residential community approximately 4 km downstream, thereby drawing significant attention [30].

2.2. Materials

The SAR data utilized in this paper were sourced from the ESA Sentinel-1A/B satellites’ images, which satellites carry C-band synthetic aperture radar. The two satellites can continuously provide SAR satellite image data for the entire world with a time interval of 6–12 days, provided that external conditions are not affected. The satellite images used in this paper are from 33 paths, with 103 views taken between 20 September 2015, and 28 December 2019, please see Table S1 for the specific data of each image acquisition. The main parameters of the SAR images are listed in

Table 1. Sentinel-1 SAR image parameters.

Orbit	Beam Mode	Repeat Period (d)	Resolution (m)	Heading Angle (°)	Incidence Angle (°)
Descending	IW	12	2.3 × 13.9	193.1	32.6

. The Guobu slope is generally laid out in the northwest direction (NW285°) with an average slope of 43° [32]. This situation means the radar pulse beam of the ascending image reaches the top of the slope first, resulting in echo signals from different areas being projected into the same SAR image cell and creating a superimposed mask. Therefore, the descending orbit data are the only option for InSAR monitoring of the Guobu slope. The external digital elevation model (DEM) utilizes the 30 m resolution SRTM3 data released by NASA, which provides accurate elevation information for the study area and improves the reliability of the deformation effect. The experiment uses precision orbit determination (POD) precise orbit ephemerides data provided by the European Space Agency for orbit refinement, effectively avoiding the systematic errors caused by incorrect orbit information.

During subsequent analysis and argumentation, precipitation and reservoir water level data were employed. The precipitation data were sourced from the meteorological station with identification number 52868, located at 36.2° N and 101.37° E, while the reservoir water level data were obtained from the Laxiwa Hydroelectric Station.

3. Methods and Materials Processing

The main research aim of this paper is to use SBAS-InSAR technology to obtain the temporal deformation sequence of active landslides and to analyze the response of said sequence to external influencing factors using wavelet tools. Figure 2 illustrates the specific workflow of this research.

3.1. SBAS-InSAR

SBAS-InSAR is a classic time series InSAR technique proposed by Berardino P et al. [33] in 2002, and it has since been further developed and refined by researchers, including Casu, F [34]. This method has been widely applied in the assessment of geological hazards such as landslides [35], earthquakes [36], and ground subsidence [37], owing to its ability to overcome the measurement errors caused by atmospheric effects, temporal and spatial decorrelation in the D-InSAR technique, and the low resolution of external DEM. The basic principle of the SBAS-InSAR method is to generate an interferogram with higher coherence by setting short spatial and temporal baseline thresholds and then performing spatial filtering (multi-looking) on the differential interferometric phases to achieve phase noise reduction. The region’s slowly-decorrelating filtered phase pixels (SDFPs), i.e., highly coherent ground target points, are identified based on the average coherence [38]. Then the observation equation is established on the identified SDFPs, and the least-squares solution of the surface deformation sequence is obtained in the sense of the minimum norm using the singular value decomposition (SVD) algorithm [33,34,39]. Finally, the atmospheric and noise phases are separated by spatial adaptive and temporal filters to obtain information such as numerical elevation model errors and temporal deformation sequences. The permanent scatterers used as part of the PS-InSAR technology are primarily located in urban areas, while SDFPs are widely distributed in natural environments, such as on exposed rocks, grasslands, bushes, etc. Therefore, SBAS is more suitable for non-urban areas such as mountainous regions and valleys [40]. Discussion on the principles of SBAS-InSAR technology has been established in numerous articles and will not be elaborated further in this paper. For further information, please refer to the citations [33], [36], and [40] in the reference list.

3.2. Wavelet Transform

The wavelet transform (WT) is a new transform analysis method that uses and develops the idea of the localization of the short-time Fourier transform, while overcoming the shortcomings of the window size that does not vary with frequency, and it can provide a “time–frequency” window that changes with frequency. This is an ideal tool for time–frequency analysis and processing. In this paper, we mainly use continuous wavelet transform, cross-wavelet transform and wavelet coherence.

3.2.1. Continuous Wavelet Transform

The continuous wavelet transform (CWT) method extends a time series into an image of spectral power in the time–frequency domain, where the wavelet power is a function of time (or date) and the period (or frequency) of the time series [41]. The resulting wavelet power image can be considered as a time–frequency map showing the signal’s power with time and frequency. Using this image, we can identify dates and corresponding periodicities that exhibit high wavelet power. These high-power regions in the time–frequency map indicate the presence of significant periodicities in the time series. Thus, the CWT method provides an intuitive way to identify significant periodicities in time series data.

Assume a discrete time series $X_n(n=0,1,\cdots ,N-1)$, which is sampled at a fixed interval $\Delta t$; its CWT can be defined as follows [26]:

$${\mathrm{W}}_i{}^X(s)=\sqrt{\frac{\Delta t}{s}}\sum _{i=0}^{N-1}{x}_i{\phi }_0^{\ast }\left[(i-n)\frac{\Delta t}{s}\right]$$

(1)

where $\ast$ denotes the complex conjugate; ${\phi }_0$ represents the mother wavelet, and in this paper, we have utilized the Morlet wavelet as the mother wavelet due to its excellent time-frequency localization capability [42]; and $s$ is the wavelet scale. The continuous wavelet transform (CWT) can be thought of as a bandpass filter with different positions and widths [43,44,45].

3.2.2. Cross-Wavelet Transform

Cross-wavelet transform (XWT) can analyze the energy resonance and covariance distribution patterns of two time series in time–frequency space, which can reveal the consistency and correlation of frequency periods of two series at different time scales [25,46,47]. Assuming two discrete time series $X_n$ and $Y_n$, the XWT between them is defined as

$$W_n^{XY}(a)=W_n^X(a)\times W_n^{Y^{\ast }}(a)$$

(2)

where $W_n{}^X(a)$ denotes the CWT of the time series $X_n$ over the frequency range $a$; $W_n{}^{Y^{\ast }}(a)$ denotes the complex conjugate of the CWT of the time series $Y_n$, and the cross-wavelet power spectrum is $\left|W_n{}^{XY}\right|$ [45,48].

3.2.3. Wavelet Coherence

As regards consistency, XWT may produce misleading results due to the fact that when the CWT spectrum of one of the time series exhibits a strong peak, the absolute value of XWT between the two will be high, even if there is no strong correlation between the two-time series. To overcome this problem, Grinsted et al. [45] proposed a normalized metric for the relationship between $X_n$ and $Y_n$, namely, wavelet coherence (WTC). WTC is defined by the following equation:

$$R_n{}^{XY}=\left|\frac{\varsigma [W_n{}^{XY}]}{\sqrt{\varsigma [W_n{}^X]\varsigma [W_n{}^Y]}}\right|$$

(3)

where $W_n{}^{XY}$ denotes the corresponding XWT, $\varsigma$ denotes the smoothing factor, and $R_n{}^{XY}$ denotes the WTC of the time series $X_n$ and $Y_n$. The value range of wavelet coherence is [0, 1], and the closer the value of WTC is to 1, the stronger the correlation between time series $X_n$ and $Y_n$ is.

3.2.4. Phase Difference and Average Phase Angle

As a complement to wavelet analysis, phase spectrum analysis can be employed to characterize correlations between signals [49]. The phase difference between time series $X_n$ and $Y_n$ can be expressed as ${\varphi }_{XY}$, reflecting the lag and relationship between the two series, enabling the determination of whether the correlation between the two is positive or negative [50]. Since the Morlet wavelet is a complex wavelet, it can be decomposed into real and imaginary components. The cross-wavelet phase difference can be defined as

$${\varphi }_{XY}=art\left[\frac{Imag\left(W_{XY}\right)}{Real\left(W_{XY}\right)}\right]$$

(4)

with ${\varphi }_{\mathrm{XY}}\in [-\pi ,\pi ]$, where $Imag$ and $Real$ denote the real and imaginary parts of the smoothed power spectrum, respectively.

We consider the eight phase angle conditions denoted by the arrows in Figure 3 [51]. If the phase difference is zero, the arrow points to the right, indicating that the changes in the time series $X_n$ and $Y_n$ are synchronized and there is a positive correlation between the two at a given frequency. On the other hand, if the phase difference is $\pi$, the arrow points to the left, indicating a negative correlation between the two. When the phase difference is $\pi /2$ or $-\pi /2$, the arrow points up or down, indicating that $X_n$ leads or lags $Y_n$, respectively. Specifically, when ${\varphi }_{XY}\in (0,\pi /2)$ and $(-\pi ,-\pi /2)$, $X_n$ leads $Y_n$, and when ${\varphi }_{XY}\in (-\pi /2,0)$ and $(\pi /2,\pi )$, $Y_n$ leads $X_n$ [52].

We quantify the phase relationship using the circular mean method, which utilizes data from regions with statistical significance above 5% beyond the cone of influence (COI). This is a valid and commonly used method to calculate the average phase of a set of angles $(a_i,i=1,\cdots ,n)$ of circular averages, and is defined as:

$$a_m=\arg (X,Y)$$

(5)

with $X={\displaystyle {\sum }_{i=1}^n\cos (a_i)}$ and $Y={\displaystyle {\sum }_{i=1}^n\sin (a_i)}$. The circular standard deviation of the wavelet phase angle is

$$S=\sqrt{-2\ln (\sqrt{X^2+Y^2}/\mathrm{n})}$$

(6)

The temporal delay or time lag (Δt) can be calculated as [41]

$$\Delta t=\frac{a_m\times T}{2\pi }$$

(7)

where $a_m$ is the mean phase angle in radians and $T$ is the periodicity or wave period of interest. Finally, the correlation of the relevant influencing factors was analyzed after calculating the hysteresis of deformation sequences detected by InSAR.

3.3. Materials Processing

In this study, Sentinel-1A/B single look complex (SLC) SAR images were processed using the standard SBAS-InSAR processing flow in SARscape software (version number of 5.6) to obtain the average deformation rate and temporal deformation sequence results of the study area. Subsequently, a deformation thematic map of the study area was produced and analyzed using ArcGIS software (version number of 10.8). This comprised the following steps:

• Preprocessing of SAR image data.

To begin, we obtain SAR images and DEM that cover the study area, and we crop them based on the research scope to shorten the processing time. In addition, the satellite orbit information needs to be corrected using the POD data, and external DEM is used for registration to improve accuracy. Before this step, the DEM in the geographic coordinate system needs to be transformed into the DEM in the SAR image projection coordinate system;

• Generation of the connection diagram.

Baseline estimation is conducted on the pre-processed SAR image data, followed by the establishment of temporal baseline connection maps using appropriate temporal baseline thresholds. For this study, the spatial baseline threshold is set at 5% of the critical baseline, and the temporal baseline threshold is set at 60 days, resulting in the generation of 419 interferometric image pairs;

• Interference processing.

Interference processing of SAR image data to obtain multiple interferograms. In this study, a multi-look ratio of 4:1 was set for the range and azimuth directions to suppress speckle noise in the communication, and the Goldstein method was employed for residual phase filtering of the interferogram. A coherence threshold of 0.2 was set, and the Minimum Cost Flow (MCF) method was used for phase unwrapping. The obtained differential interferograms and unwrapped phase maps were visually interpreted, and 202 pairs with unsatisfactory interferometric and unwrapping effects were eliminated. Finally, the spatio-temporal baseline connection map is presented in Figure 4;

• Refinement and reflattening.

Orbit refinement and phase shift calculation are performed using cubic polynomials based on the selected ground control points (GCPs) to eliminate possible slope phases and correct for satellite orbits and phase shifts;

• Deformation inversion estimation.

The first inversion is conducted to estimate the deformation rate and residual topography, and the corresponding wavelet decomposition level is set to remove the residual topography as needed. Subsequently, the second inversion is performed, with high-pass filtering in the time domain and low-pass filtering in the spatial domain, to eliminate the influence of atmospheric phase. Finally, the time series deformation results are obtained;

• Geocoding.

The deformation results in the SAR coordinate system are converted to the geographic coordinate system of the reference DEM, and then the Arcgis software is used to create thematic maps, such as annual average surface deformation rate and cumulative deformation maps, for the subsequent related analysis.

The time series used for wavelet analysis requires equally spaced temporal sampling. Despite the fixed 12-day revisit period of the Sentinel-1A/B satellite, some images may be missing or excluded from the SBAS-InSAR processing stream, resulting in a nonhomogeneous temporal distribution of the LOS deformation time series. Therefore, we used linear interpolation to fill in missing deformation data. It should be noted that linear interpolation does not introduce false periodicity, as its linear features only produce a constant or DC component in the frequency domain [47]. The rainfall and water level data collected in this study were sampled daily. To maintain the same sampling interval as in the InSAR time series, we downsampled the rainfall and water level data. However, due to the differences in the physical meanings between rainfall and water level data, we treated the cumulative rainfall over the previous 12 days as the rainfall data for each point when downsampling the rainfall data, and the average water level over the previous 12 days as the water level data for each point when downsampling the water level data.

4. Results and Discussion

4.1. SBAS-InSAR Results

To investigate the surface motion status of the Guobu slope during 2016–2019, this study detected the surface deformation of the Guobu slope using the SBAS-InSAR technique. Based on the method described in Section 3, a time series deformation sequence of Sentinel-1A/B SAR image data was generated, resulting in a total of 28,498 deformed pixels in the studied region. The relationship between deformation rate, deformation rate precision, and pixel number is illustrated via the histograms in Figure 5. Figure 5a shows that the histogram is negatively skewed with a mean deformation of −45 mm/y, indicating that some regions in the studied area are subsiding. More than 75% of the monitored points have deformation rates distributed in the range −20 to 20 mm/y, indicating that most of the studied region is stable. Figure 5b shows that the precision of the pixel deformation rate in the entire study area is in the range of 1 to 8 mm/y, with a standard deviation of 0.967 mm/y. These two histograms are, to some extent, indicative of the validity and accuracy of the experimental results.

We utilized ArcGIS to produce an annual average deformation velocity map of the Guobu slope and present it in Figure 6. Figure 6a exhibits the annual average deformation velocity from September 2015 to December 2019, while Figure 6b displays the cumulative deformation map at certain moments. Negative values indicate that the landslide body is moving away from the SAR sensor and vice versa. We calculated the annual cumulative deformation at the largest deformation point and listed the results in

Table 2. The annual cumulative deformation during the monitoring period.

Time Period	Cumulative Deformation/mm
2016.1~2016.12	−430
2017.1~2017.12	−445
2018.1~2018.12	−442
2019.1~2019.12	−433

. As illustrated in Figure 6, the maximum deformation rate of the Guobu slope during September 2015 to December 2019 reached −447.448 mm/y, and the maximum cumulative deformation was −1860 mm. The landslide boundary was distinct and shows a “U” shape in the plan view. These results are consistent with the monitoring results of Shi et al. [30] and Li et al. [27] in terms of spatial distribution and deformation magnitude, further confirming the reliability of our deformation monitoring results. Table 2 indicates that the movement of the Guobu slope is stable at the interannual scale, without apparent acceleration or deceleration (the 15 mm annual deformation difference can be ignored compared to the 440 mm cumulative deformation).

In order to comprehensively investigate the spatial deformation characteristics of the Guobu slope, a profile line denoted as AA′ is drawn from the crest to the toe of the slope along the primary sliding direction of the deformed mass (depicted as a solid black line in Figure 6a). Subsequently, the point data collected along the profile line AA′ are depicted in the form of a polyline graph, exhibiting the relationship between the deformation rate and elevation in the direction of the profile line, as illustrated in Figure 7. The solid blue line in Figure 7 represents the elevation along the profile line AA′, whereas the solid red line represents the annual average deformation rate along the profile line. The dashed green line corresponds to the elevation of the design water level at 2452 m, and the dashed black line serves as an auxiliary guide to determine the elevation of the severely deformed region. Notably, the region of most severe deformation, where the annual mean deformation rate exceeds −400 mm/a, is located in the upper middle section of the slope, between 2700 and 2850 m ASL. In addition, the entire slope has an annual mean deformation rate of more than −200 mm/a. The annual deformation rate from the top to the bottom of the landslide body demonstrates a “V”-shaped variation, transitioning from gradual to quick to gradual, while the displacement of the slope surface exhibits a typical nonlinear spatial and temporal distribution.

We selected three distinct target points on the slope surface for time series analysis, as shown by the blue dots in Figure 6a. To mitigate the interference of random errors on these selected points, a 60 m buffer zone was created around each target point. The mean displacement of all deformation points within the buffer zone was chosen as the representative value for each target point, as shown in Figure 8. The error bars indicate the standard deviation of all measurements within the buffer zone. Notably, point P1, located in the upper-middle part of the landslide, has an accumulated deformation of over 1.75 m during the monitoring period, while point P2, situated in the bottom area of the landslide, has an accumulated deformation of approximately 1.2 m. Point P3, located in the central area of the plateau, has an accumulated deformation close to 1 m.

The deformation time series of the feature points exhibits non-uniform linear changes from September 2015 to December 2019, which can be regarded as the arithmetic sum of linear and nonlinear components. The linear component represents the non-elastic displacement of the slope due to its own gravity, while the nonlinear component represents the acceleration or deceleration of displacement due to external factors [53]. In this study, the linear component of deformation was obtained through least-squares fitting, and the nonlinear component is defined as the difference between the deformation time series and the linear component. A plot of the displacement nonlinear component versus precipitation and water level data is constructed and shown in Figure 9. It can be observed that there is a certain periodicity and correlation between the nonlinear component of the feature points and the precipitation and reservoir water level data. To additionally discover their correlations, the wavelet transform method is employed to analyze their correlations and lags in the time–frequency domain.

We extracted the reservoir water level data from 2009 to 2014 and the 3D deformation data of 17 GPS monitoring points located in the headscarp of the Guobu slope from Zhang Tao’s article [26], and used vector displacement, i.e., $\sqrt{N^2+E^2+U^2}$, to represent the landslide movement, as shown in Figure 10. Shi et al. performed a cross-correlation analysis of the reservoir water levels and the nonlinear components of 15 GPS monitoring points (excluding QC9 and IP10) between June 2011 and November 2013. They found that the changes were fully synchronized, indicating that the stability of the slope was entirely influenced by the fluctuation of the reservoir water levels and that the impact of precipitation could be disregarded [30].

4.2. Wavelet Transform Sesults

4.2.1. CWT Result

We conducted CWT on the rainfall, reservoir water level, and nonlinear deformation components of feature points to show their time–frequency features. The CWT is expressed as a two-dimensional graph, i.e., continuous wavelet power spectrums with two axes as the proxies of time instant and period (or frequency) of the specific pattern (Figure 11). Therefore, a location with a high modulus value in the CWT generally indicates a localized signal that significantly resembles the frequency-specific wavelet.

It has been ascertained from the continuous wavelet power spectrum that there is a period of significance of approximately 1 year (30 × 12 days) for both rainfall and feature point P2 throughout the temporal domain of 2016 to 2019, as depicted in Figure 11a,d. Furthermore, the reservoir water level demonstrates a period of significance of approximately 1 year between January 2017 to August 2019, as illustrated in Figure 11b. Similarly, feature point P1 displays a prominent period of around 1 year from January 2016 to July 2017, as portrayed in Figure 11c. Compared to other data, it should be noted that feature point P3′s periodicity of around 1 year is not particularly pronounced, even though there is noticeably high energy density in the power spectrogram (Figure 11c).

In Figure 11a, the peak signal exhibiting a periodicity of 0–48 days appeared in June of 2016, 2017, and 2018, as well as in September of 2019. This may be related to the seasonal precipitation occurring predominantly during summer in the Laxiwa region. A similar phenomenon appears in Figure 11b, where the peak signals of days 0–48 and 72–120 occurred around June of each year, which are likely related to the fluctuations in the water level of the reservoir that could themselves have been caused by flood control, power generation, and irrigation operations at the Laxiwa hydroelectric station. Furthermore, notable power signals with a periodicity of 1–3 months are observed in Figure 11c–e. The underlying causes of these signals require further explanation.

4.2.2. XWT and WTC Results

The XWT and WTC techniques exhibit strong signal coupling and resolution capabilities, enabling them to identify common power and relative phase information in the time–frequency domain and assess the significant coherence and confidence levels against the red noise background [25]. By performing XWT and WTC transformations on the two time series obtained from CWT, we obtain the cross-wavelet energy spectrum and the cross-wavelet condensation spectrum between the nonlinear deformation components of the characteristic points and the rainfall and reservoir water levels (Figure 12 and Figure 13).

In Figure 12a,c,e, we can observe that the cross-wavelet energy spectrum between precipitation and three feature points exhibit the highest energy density in a 1-year (30 × 12 days) period and passed the 95% significance test, with significant periods continuously distributed in the entire sequence domain. The cross-wavelet arrows between precipitation and feature points P1 and P2 are directed toward the lower right (Figure 12a,c), indicating that the rainfall is ahead of the nonlinear displacements. This fact indicates that the nonlinear displacements exhibit a time lag response to the rainfall. To quantify the time lag within the interannual structure, we calculate the average phase angle and circular standard deviation from the phase relationship diagram of Figure 3 and Equations (5) and (6), and calculate the corresponding time delay days using Equation (7) (Shown in

Table 3. Average phase angle and time lag of rainfall, reservoir water level, and nonlinear displacement.

Factors	Point	Significant Period/d	Significant Time/year	Average Phase Angle/rad	Time Delay/d
Rainfall	P1	314.8~396.7	2016~2019	−1.183 ± 0.210	68.7 ± 12.2
	P2	297.8~420.3	2016~2019	−1.246 ± 0.152	72.4 ± 8.8
	P3	405.6~428.3	2016~2019	1.433 ± 1.887	83.3 ± 109.6
Reservoir water level	P1	353.4~396.7	2016~2019	1.762 ± 0.313	80.1 ± 18.2
	P2	280.5~420.3	2016~2019	2.810 ± 0.142	19.3 ± 8.3
	P3	—	—	—	—

). The deformations of P1 and P2 lag behind precipitation by 68.7 ± 12.2 and 72.4 ± 8.8 days, respectively. This finding has significant implications for our understanding of the complex interplay between precipitation and deformation in this particular area. The cross-wavelet arrow between precipitation and P3 is unstable, with a circular standard deviation greater than the average phase angle, indicating that the two do not have a stable correlation, probably due to the impact of human factors on P3, located on the headscarp of the Guobu slope (Figure 6A), which require further investigation.

In Figure 12b,d, we can observe that the cross-wavelet energy spectrum between reservoir water level and feature points P1 and P2 exhibit the highest energy density in a 1-year (30 × 12 days) period, and pass the 95% significance test, with significant periods continuously distributed in the entire sequence domain. The arrows point upward to the left, indicating a negative correlation between P1 and P2 and the reservoir water level. In other words, when the reservoir water level decreases, the deformation of P1 and P2 increases. According to Table 3, the deformations of P1 and P2 lag behind the water level by 80.1 ± 18.2 days and 19.3 ± 8.3 days, respectively. The time lag at P2 is significantly smaller than that at P1 because P2 is located at the foot of the slope (Figure 6A) and is closer to the reservoir water level. However, Figure 12f shows no effective resonance period between the reservoir water level and P3.

The WTC results reveal correlation and time lag information that is consistent with the XWT results, although there are some differences in local details. Figure 13a,b show the WTC between rainfall and nonlinear displacement P1 passes the 95% significance test in the 1-year (30 × 12 days) period, while the WTC between reservoir level and P1 do not. This result implies that the nonlinear displacement of P1 is mainly affected by precipitation from 2016 to 2019. Figure 13c,d show that the WTC results for rainfall, reservoir water levels, and P2 pass the 95% significance test in the 1-year (30 × 12 days) period. However, the color of the wavelet coherence spectrum indicates that the WTC value between precipitation and P1 (about 0.85) is significantly higher than that between reservoir water level and P1 (about 0.7), suggesting that the nonlinear displacement of point P2 is influenced by both precipitation and reservoir water level at a period of 1 year, with rainfall having a greater contribution. In Figure 13e,f, there is no apparent coherence between P3 and precipitation or reservoir water level in the 1-year period. Still, a coherence patch is found between the two images at the intersection of the horizontal axis of 2018 and the vertical axis of the 45–120-day period, which may be related to the large fluctuations in precipitation and reservoir water level in 2018.

4.3. Challenges While Processing

This article and the related research [12,13,36,39,53] have proved that the time-series InSAR technique is an effective and efficient tool for landslide monitoring, but it suffer from certain limitations and shortcomings. Landslides are typically found in complex topographic regions, such as mountains, hills, and river valleys, where SAR side-looking imaging mode can result in SAR images having overlapping, shadowing, and foreshortening. When employing time-series InSAR for landslide monitoring, some landslides may not be detected. Additionally, the sliding direction of landslides can also cause a decrease in the monitoring accuracy of InSAR technology, particularly when the direction of landslide movement is perpendicular to the satellite’s LOS, causing the displacement component of the landslide to be zero in the satellite’s LOS, thereby rendering it impossible for time-series InSAR to detect the landslide movement. Therefore, when employing the time-series InSAR technique for landslide monitoring, it is necessary to take into account various factors such as terrain, landslide motion direction, and SAR satellite line of sight direction. In exceptional circumstances, a combination of satellite data from both ascending and descending track or multiple satellite platforms can be utilized.

Similarly, we encountered some issues when performing wavelet transforms on our time series data: (a) wavelet tools result in an edge effect at the beginning and end of the time-frequency space (the cone of influence), which cannot be eliminated and decreases the spectral accuracy of landslide time-series information [47]; (b) XWT and WTC require equal time intervals for the time series, which reduces the temporal resolution of certain time series; (c) XWT and WTC can only analyze two time series, reducing the computational efficiency when multiple time series are involved, and making it difficult to explain the correlation and time lag information between multiple series. In the future, we could explore improved methods, such as least-squares wavelet (LSWAVE) [23], multiple wavelet coherence (MWC) [54], etc., or other methods in time series analysis [55], to solve the above problems of InSAR time series with drivers analysis.

5. Conclusions

Utilizing the SBAS-InSAR technology, this study closely monitored the active slope upstream of the Laxiwa Hydropower Station, yielding valuable insights into the displacement and deformation time series of the Guobu slope from 2016 to 2019. Notably, the resulting high-precision landslide displacement data proved effective in recording the movement with accuracy and reliability that were verified using self-data accuracy assessment and previous monitoring results. The InSAR monitoring results of the time series reveal that the steeper area in the middle and upper parts of the Guobu slope display the largest deformation rate, validating the failure mechanism of the tilting deformation of the Guobu slope. By incorporating precipitation, reservoir water level data, and the deformation sequence obtained through InSAR technology, this study quantitatively analyzed the response of landslide movement to hydrological changes using wavelet tools. The findings indicate that seasonal precipitation predominantly influenced landslide movement during 2016–2019, and the deformation responded to precipitation with a delay of approximately 70 days. Additionally, the lower part of the slope exhibited a quicker response to water level changes at about 20 days, while the area with the largest deformation in the middle and upper parts of the slope responded to water level changes at around 80 days. The landslide movement is currently mainly driven by precipitation, compared to the 2009–2014 reservoir storage period during which the landslide movement was mainly driven by the reservoir water level. Through the application of wavelet tools, this study quantitatively explains the correlation and lag between landslide movement and influencing factors. This study thus provides a new approach to the assessment of landslide movement and facilitates future disaster prevention and control.