Stability analysis of reference station and compensation for monitoring stations in GNSS landslide monitoring

The Real-Time Kinematic (RTK) positioning method of the Global Navigation Satellite System (GNSS) has been widely used for landslide monitoring. The stability of its reference station is crucial to obtain accurate and reliable monitoring results. Unstable reference stations due to the geological environment and human activities are difficult to detect and in practical applications often ignored. As a result, it affects the positioning solutions and subsequently the interpretation and detection of landslide motions, which must be addressed in GNSS landslide monitoring. To solve this problem, we propose using the Precise Point Positioning (PPP) technique to analyze the stability of the reference station by verifying its position. The deformations of the monitoring stations are then compensated. First, the reference station coordinates are obtained by the PPP technique and tectonic motion is considered in data processing. The change or breakout of the reference station position is then determined using a cumulative sum control chart method. Finally, each monitoring station’s displacements are compensated according to the displacements of the reference station. According to the results of the Tengqing landslide experiment, the PPP technique can be used in GNSS landslide monitoring to analyze the stability of reference stations. With PPP, millimeter-level accuracy for the coordinates of reference stations is achieved. Compared to the traditional deformation series, the compensated displacement series more reliably reflects the landslide motions. This study will increase the reliability of monitoring results and contribute to implementing GNSS in monitoring landslides.


Introduction
Many countries have been facing various geological disasters, particularly landslides.(Froude & Petley, 2018).Landslides cause a great loss of live and property every year (Cruden & Varnes, 1996;Yaghmaei, 2020).Its movements must be continuously monitored for effective disaster prevention (Tomás & Li, 2017).The Global Navigation Satellite System (GNSS), which has advanced quickly in recent years, has emerged as one of the most accurate techniques to monitor landslides (Benoit et al., 2015;Jing et al., 2022;Tu et al., 2017;Zhang et al., 2022).For instance, the Variometric Approach for Displacement Analysis Stand-alone Engine (VADASE), developed by Leica, is a processing approach capable of estimating the velocities and displacements in real-time using GNSS data in a global reference frame (Benedetti et al., 2014).RT-SHAKE, an algorithm presented by Psimoulis et al. (2018), is designed to detect ground motions related to various natural phenomena.In addition, Bai et al. (2019) developed a real-time BeiDou Navigation Satellite System (BDS) monitoring and warning cloud platform.This platform successfully applied early warning techniques to the Heifangtai landslide in Gansu Province, China.
Absolute positioning and relative positioning are the two primary techniques for processing GNSS data (Wang, 2013).The Real-Time Kinematic (RTK) technique, which uses the accurate measurements collected simultaneously at a reference station and a monitoring station, is the method used to implement the relative positioning strategy.Based on this, a relative position with respect to a reference station can be estimated (Du et al., 2019).This method can get rid of orbital error, tropospheric, and ionospheric delays when the baseline is short (Han et al., 2018).The absolute positioning approach such as Precise Point Positioning (PPP) can be completed with a receiver (Huang et al., 2017).While PPP can achieve positioning with centimeter-level accuracy but needs a lengthy convergence period of about 30 min (Li et al., 2015).
The RTK technique has been widely used to precisely monitor landslide motions.This technique typically installs a stable reference station outside the deformation area and provides the correct information for one or more monitoring stations (Wang & Soler, 2012).The position coordinates of a reference station, which remain unchanged over time, are used to compute relative displacements of the monitoring stations.However, there may also be the external influences that can inadvertently change the position of the reference station, such as plate movement, subsidence, creep area spread, and human activities in urban environments.These influences will have a detrimental effect on the reliability and continuity of the positioning results for the reference station (Bao et al., 2018;Wang et al., 2019Wang et al., , 2022)).Such external influences, which are difficult to detect and avoid and resolve, must be taken into account in GNSS landslide monitoring.
One solution is to use PPP technology for absolute positioning (Lin et al., 2021;Martín et al., 2015;Yigit et al., 2016).Capilla et al. (2016) studied the PPP technique for landslide monitoring using Global Positioning System (GPS) and GLObal Navigation Satellite System (GLONASS) combined.This study concludes that the real-time PPP technique can measure deformations with a standard deviation of only 2 cm, which ensures that the real-time PPP technique is capable of detecting deformations up to the magnitude of 5 cm.Wang et al. (2013) defined a stable reference frame using the data from 8 permanent stations in Puerto Rico and the Virgin Islands, and realized millimeter-level accuracy of monitoring stations based on PPP technology.However, limited by the accuracy and real-time performance of PPP technique, this method is typically used for monitoring low-rate, daily geological hazards for a long period (Huang et al., 2023).
Another solution is to use the RTK technique to determine the three-dimentional displacements.Additionally, the stability of the reference station is monitored using absolute positioning.Du et al. (2020) proposed and applied an asynchronous RTK method to determine if there is a measured shift at the reference station and provide accurate movement results.However, the motion of the reference station is typically long-term, continuous, and at low-rate, making it only for detecting the displacements that occur over a short period and not suitable for long-term stability monitoring of the reference station.Bian et al. (2014) established a high-precision GNSS monitoring network using the International GNSS Service (IGS) stations in the vicinity as reference points.This method can be used to obtain the reference station coordinates.However, this method is also dependent on the reference stations, which are restricted by the length of the baseline.Data processing based on the PPP technique represents a possible solution to obtain reference station coordinates without any direct cross-correlation between stations.This technique has been well-established and frequently used in the field of ground observation, allowing us to derive reliable and consistent results.Although the PPP technique cannot get high-accuracy coordinates in real time, but due to the daily displacement of the reference station is small and negligible, the daily and post PPP technique is sufficient to meet the stability analysis of the reference station.
Aiming at the long-term, continuous, and slow movement of a reference station, this paper proposes to use the PPP technique and control chart method to analyze the stability of the reference station and to compensate for the monitoring station's displacement.An experiment on the Tengqing landslide demonstrated the feasibility and practicality of the proposed technique, assuring early warning timeliness and monitoring accuracy.The remaining sections of the paper are organized as follows.Sect."Methodology" introduces the stability analysis and compensation methods for the reference station.Sect."Experiment and results analysis" presents the experimental results on the Tengqing landslide to support the proposed strategy.The final section of this study presents our conclusions.

Acquisition of reference station coordinates
Obtaining accurate coordinates of a reference station is a prerequisite for the stability analysis.PPP is used for the non-differenced observations when the GNSS receivers observe the pseudo-ranges and carrier phases.(Wang et al., 2021).Moreover, the IGS's precise satellite orbit and clock products can be used to obtain the best position solution (Dow et al., 2009;Héroux & Kouba, 2001;Xie et al., 2023).GNSS derived position is in an Earthcentered, Earth-fixed three-dimensional Cartesian coordinate system (ECEF-XYZ).In static mode, the PPP approach delivers the position with millimeter to centimeter precision, which is used to obtain the ECEF-XYZ coordinates (X(t) IGS14 , Y (t) IGS14 , Z(t) IGS14 ) at epoch t in the source reference frame IGS14.
Regarding the RTK technique, the reference station is close to the study area and exhibits minimal tectonic and geophysical loading effects.This ensures that the local landslide movements can be captured more accurately.This is a significant advantage of the RTK technique.However, the PPP technique may be influenced by tectonic motion and geophysical loading, leading to potential bias in the stability analysis of the results.A local reference frame can be used to overcome this problem.In general, Helmert transformation is used to implement local reference frames.Wang et al. (2018) used a seven-parameter transformation approach to convert the coordinates from the IGS global reference frame to the stable North China Reference Frame 2016 (NChina16).The equations for obtaining ECEF-XYZ coordinates in a regional reference frame are as follows: where X(t) R , Y (t) R , and Z(t) R are the station's XYZ coor- dinates in the target reference frame at epoch t , t 0 indi- cates a specific epoch that is utilized to link the source reference frame (IGS14) with the target reference frame, and the abbreviations T x ′ , T y ′ , T z ′ , R x ′ , R y ′ , and R z ′ stand for the six transformation parameters.
Seasonal variations are not included in the local reference frame.Hydrology and atmospheric loading are frequently regarded as the main causes of the variations in these areas.This frame realization may be biased by seasonal variations in the frame sites, which are aliased into the site coordinates.Since this seasonal signal might affect the characteristics of interest that are inferred from (1) the time series, it is crucial to model it (Davis et al., 2012).Interseasonal variability has been considered in the earlier analyses.Freymueller (2009) assessed seasonal fluctuations at a number of sites in northern North America and the Arctic using an iterative non-parametric technique.At most locations, seasonal changes in the vertical surpass 4-5 mm.Zou et al. (2014) constructed seasonal loading models for each site from the time series of loading displacements projected by either the Gravity Recovery and Climate Experiment (GRACE) load model or the collection of forward models.Additionally, the alignment of the GPS solution to an enhanced Terrestrial Reference Frame solution that takes the loading model into account was used to account for the seasonal fluctuations in loading.Consequently, a linear + seasonal model can be considered as a component of the reference frame.

Stability analysis of reference station
After obtaining the coordinates of the reference station using the PPP technique and mitigating the tectonic motion using the local reference frame, non-tectonic and non-landslide deformations are left.These deformations can result from various natural and anthropogenic (human-induced) processes and are the main cause for reference station instability.Examples include glacial isostatic adjustment, post-glacial rebound, sediment loading and compaction, groundwater extraction, volcanic deformation, reservoir-induced deformation, and human activities.The distinction between these non-tectonic and non-landslide deformations depends on the specific geological, geographical, and anthropogenic factors in a given area.It is necessary to analyze the effects of nontectonic and non-landslide deformations on the stability of the reference station and to mitigate them.
Displacement in the horizontal and vertical directions is of greater concern.Thus it is necessary to convert the coordinates as the origin of the coordinate system. (2) where B(t 0 ) is the longitude of the initial coordinate and L(t 0 ) is the latitude.For the convenience of analysis, the three displacement components are converted into one quantity as Eq. ( 3) because it represents the landslide displacement more significantly. (3) The next step is to determine the changes or breakouts at the reference station by CUmulative SUMs (CUSUM).An approach to achieve this is a CUSUM control chart method (Hawkins & Olwell, 2012).The first 25 samples of the sequence x 1 , x 2 , x 3 , . . ., x n are used to estimate the target average m x and target standard deviation σ x .Two cumulative sums are tracked by CUSUM: The upper cumulative sum detects when the local mean shifts are upward: The lower sum detects when the mean shifts are downward: The variable s is the minimum detectable mean shift with a default number of 1.If a process follows U j > cσ x or L j < −cσ x , it violates the CUSUM condition at sample x i , where the default value for the variable c is 5, which indicates how many times of the standard deviations from the mean the upper and lower cumulative sums are allowed to depart (Chang & McLean, 2006;Sibanda & Sibanda, 2007).As a result, CUSUM examines the data for the first 25 samples and alarms whenever there is a mean shift of more than five standard deviations from the initial data.

Influence of reference station displacement on monitoring station
The undifferenced carrier phase observation for the satellite i and receiver A can be expressed as (Liang et al., 2015): (4) where the abbreviations ( ϕ , ρ , , N , c , I , T , E , ε ) repre- sent the carrier phase measurement, geometric distance, signal wavelength, carrier phase ambiguity in cycles, light velocity, ionosphere delay, troposphere delay, and measurement noise of the carrier phase and multipath.For a short baseline, the difference equation between the monitoring receiver A and reference receiver B elimi- nates or attenuates some errors (Takasu et al., 2010): is the direction cosine between the station's approximate location and satellite i.
Assuming the reference station B is displaced by m along the X direction.Since the baseline length has not changed, the X-direction component ( l i B dX B − l i A dX A ) in Eq. ( 7) remains unchanged.To ensure the establishment of Eq. ( 7), the monitoring station should also be displaced by n: where n is The value of n depends on m and the direction cosines at the baseline's two ends.The longer the baseline length, the greater the difference between l i B and l i A will be.In the conventional landslide monitoring scene, the baseline length is generally short (< 10 km), so the difference between l i B and l i A can be ignored.As a result, the displacement of the monitoring station n is equal to the amount of motion m of the reference station.

Compensation for the monitoring station displacement
According to Eq. ( 9), when the reference station's coordinates remain unchanged, the calculated displacement is biased by m in the opposite direction from the calcu- lated monitoring station's displacement.In this way, the monitored displacement at each monitoring station can be compensated after the displacement at the reference station is obtained.The compensation model is: (6) where (N (t) A )′ , (E(t) A )′ , (U(t) A )′ represent the moni- tored displacement after compensation.However, the PPP technique is prone to errors due to various factors such as uncorrected residual errors, monitoring environmental changes, convergence time, atmospheric delays, and receiver noise.These factors can lead to missing or noisy data which do not accurately reflect the true deformation.Preprocessing of the landslide displacement series is a critical step to reduce the factors of these noise sources.We proposed using the K-Nearest Neighbors (KNN) algorithm to fill in missing data.The purpose of the KNN algorithm is to find 'k' samples in the dataset that are comparable or nearby in space (Lee & Styczynski, 2018).These 'k' samples are then used to predict the values of the missing data points.For each sample, the missing values are imputed using the mean values of the k-neighbors that are obtained in the dataset.The LOcally Weighted Regression (LOWESS) algorithm was used to smooth the PPP results.The samples are primarily divided into brief intervals, and the samples within the interval are fitted with polynomials (Royston, 1992).To get weighted regression curves for various periods, this procedure is repeated.In order to create a full curve, the centers of each regression curve are finally connected.In addition, other methods can also be used to enhance the reliability of the analysis, such as spatial regression to evaluate the spatial correlation of stations (Habboub et al., 2020).

Processing flowchart
The technical flowchart of stability analysis of reference station and compensation for monitoring stations in GNSS landslide monitoring is depicted in Fig. 1.Three steps make up this data processing procedure.(2) The reference station's cumulative displacement in the one-dimensional direction is obtained by the coordinate transformation, and the status of the reference station is analyzed using the CUSUM method.An alarm is raised when the mean change from the initial data reaches five times of standard deviation.
(3) Each monitoring station's displacement series is adjusted when the reference station's instability is discovered in accordance with the idea that the monitoring station is affected by the reference station's shifting.To obtain higher accurate compensation values for the reference stations, we use the KNN algorithm to compensate for missing values and the LOWESS algorithm to smooth the noisy results.

Experiment area
This experiment was conducted in the coal mining district of Liupanshui, Guizhou Province, Southwest China, in the area of the Tengqing landslide in Fig. 2. The Tengqing landslide study area is indicated by the red polygon.There are residential properties along the edge of the road at the bottom of the landslide area, with forest covering the top.The annual average temperature of the experiment area is 12.3 °C and the rainfall is 1 100-1 300 mm, which experiences a plateau monsoon climate.Since 2011, underground coal mining has significantly deformed the landslide area, and the development continues at a rapid pace.A number of stones fell to the ground after being loosened from the rock's face.The residents in the landslide area were relocated to prevent property and life loss (Ren et al., 2022).There are still some factories and workers in the vicinity of the landslide, and the mine continues to operate under normal mining conditions.To learn more about the internal workings of landslides five lowcost GNSS monitoring stations were installed in October 2019.Among them, a reference station (TQ01) was installed outside the landslide which was stable in the project design phase and had a good observation environment.Considering the possible subsidence caused by coal mining, a backup reference station (TQ02) was installed away from the monitoring area.And this station was considered as a monitoring station in this paper.Three monitoring stations were installed on the landslide deformation body (from TQ03 to TQ05).The distribution of these GNSS stations in the Tengqing landslide is shown in Fig. 2.

Data processing with a stable reference station
The evaluation of the proposed approach was conducted using GNSS monitoring data from October 4, 2019 to October 1, 2020.A GNSS software program developed by our research team was used to process the data.Our experiments show that it can monitor the landslide movements at millimeter-level accuracy in an open environment (Bai et al., 2019;Huang et al., 2018).A partial ambiguity fixing technique was used in the BDS/GPS RTK positioning solution.The broadcast ephemeris was used, and the satellite cut-off height was 10°.The parameter setting of the GNSS software and station details are given in Table 1.The distances between the reference stations and the monitoring stations are represented by the baseline length in the table.
For ease of analysis, the daily deformation characteristic was derived from the median of the real-time values (Huang et al., 2022).Figure 3 shows the calculated GNSS displacement series for the four stations from TQ02 to TQ05.According to the displacement series, all four stations experience large displacements, and the displacement trends are consistent.Moreover, the direction of the displacement is upward, which is inconsistent with the evolution of landslides.Therefore, it is crucial to assess the stability of the reference station if RTK is used for landslide monitoring.

Stability analysis of the reference station
We processed the data at the reference station using RTKLIB software in a GPS PPP strategy (Du et al., 2021;Takasu & Yasuda, 2009).We used GPS dual frequency observations with a 10-degree elevation cutoff and the IGS14 reference frame's high-precision GNSS orbit and clock products.Besides, we used the Earth orientation parameters from IGS.We used the Zenith Tropospheric Delay (ZTD) model to estimate the tropospheric delay and the ionospheric-free method.The parameter setting of the PPP technique for the reference station TQ01 is given in Table 2.We used the Stable South China Reference Frame 2020 (SChina20) established by Bao et al. (2021) to avoid the biases caused by plate motion.Table 3 lists the SChina20 implementation parameters and their units.
The displacement series of reference station TQ01 is shown in Fig. 4. The first 25 days of data without significant deformation were selected to evaluate the results of PPP.In the N, E, and U directions, the coincidence accuracy is 1.9 mm, 5.1 mm, and about 6.1 mm, respectively.According to the positioning results, the cumulative deformation in the N, E, U, and R directions reached 185.6mm, − 95.6mm, − 875.8mm, and 900.5mm, respectively.The reference station is unstable and it is necessary to determine the time when the instability occurs.If the effects of tectonic movement and geophysical loading are not considered, the cumulative deformation in the N, E, U, and R directions reached 183.6 mm, − 88.7 mm, − 875.8 mm, and 899.2 mm, which will have a deviation of several millimeters.Therefore it is better to transform the positioning results from the IGS global reference frame to the stable local reference frame.
The mean and standard deviation were calculated for the first 25 days in the R direction, which is 8.0 mm and 3.8 mm.They serve as the target mean and standard deviation in the CUSUM control chart.Figure 5 indicates the results of the stability analysis applied to  the CUSUM control chart in the R direction.The range of dates is from October 4, 2019, to December 22, 2019, and the vertical axis indicates the standard error, which is the ratio of the cumulative total to the desired standard deviation.The locations where the cumulative sum from the target mean is greater than five times of standard deviation are highlighted.It can be observed from Fig. 5 that the lower sum does not drift more than five times of standard deviation, indicating that the reference station is not deformed in the opposite direction.The drift result of the upper sum exceeds five times of standard deviation on November 28, 2019, when the reference station instability is detected.Compared to the traditional threesigma criterion, CUSUM effectively enhances the ability of the GNSS landslide monitoring systems to identify the instability of the reference stations.

Displacement compensation for each monitoring station
After detecting the instability of the reference stations, the monitoring results of each monitoring station can be compensated.As depicted in Fig. 4, the 2020-05 and 2020-06 data feature a lot of noisy and missing data.
We use the KNN algorithm to compensate for missing values and the LOWESS algorithm to smooth the noisy results.The k-value of the KNN algorithm used to fill in the missing values is 10, and the interval width of the LOWESS algorithm used for smoothing is 3% of the observed counts.Figure 6 shows the series of the reference station TQ01 after filling in the missing values and smoothing.Although there are differences between the sequences before and after preprocessing, it is important to note that the overall trend remains consistent.The sequence of reference station TQ01 after preprocessing is close to the real displacement and could be used to compensate for the displacement results of the monitoring station according to Eq. ( 9).The displacement time series of TQ01 served as the basis for the compensation terms used for displacement correction.
The series for the four monitoring stations, both before and after compensation, are also shown in Fig. 6.The detailed cumulative displacement and the average rate for each station are given in Table 4.The average deformation rate of TQ02 after compensation is less than 0.05 mm/d and there are no significant deformations, which is consistent with the actual scenario.This consistency proves that the method effectively compensates for the shift caused by the motion of the reference station.The studied area has experienced both the subsidence of the mining area and landslide evolution, according to the findings of the field exploration.
The local geological features, mining practices, and other factors can greatly influence the magnitude and direction of the deformations observed at different stations.The subsidences at TQ03, TQ04, and TQ05 were particularly significant and aligned with the estimated deformations.However, there are variations in the displacement directions among these stations in the north and east.This inconsistency can be attributed to the different impacts of the two disasters on different stations.At present, the three stations undergo a uniform deformation, but continuous tracking of their future motions is required.
Since the real displacement of the monitoring station in the first 25 days is consistent with the traditional displacement sequence, we selected that specific period to evaluate the accuracy of displacement compensation.The Root Means Square Error (RMSE) in the N, E, and U directions are 1.9 mm, 3.1 mm, and 4.6 mm, respectively, satisfying high-precision landslide monitoring.The above explanation indicates that the proposed method can be utilized to account for the monitoring station offsets caused by the reference station's instability and provide reliable deformation results, which is essential in landslide monitoring.

Evaluation with the backup reference station
TQ02 is a backup reference station located far from the monitoring area.The compensation results in the previous section prove that TQ02 is stable.Therefore, we can use the backup reference station TQ02 to verify the displacement of TQ01 and whether the compensation results are reliable.We evaluated the solution using TQ02 as the reference station and obtained the results, which are presented in Fig. 7.The displacement trends of each monitoring station when TQ02 is used as the reference are consistent with the PPP and compensated results.This illustrates that the challenges associated with high-precision monitoring when the reference station is unreliable are effectively addressed by the compensatory strategy suggested in our research.
We chose a landslide in the upstream of the Three Gorges Reservoir to further confirm the method's general applicability.The evaluation was performed using GNSS monitoring data from January 1, 2020 to October 1, 2020.
The landslide has a substantial displacements during the flood season, which is caused by rainfall during the rainy season and fluctuations in the water level of the Three Gorges Reservoir.This can also affect the reference station.Table 5 shows the results of applying the method to compensate each monitoring station (XP06-XP10), which enhances the reliability of the landslide monitoring.

Conclusions
Landslide deformation can be accurately monitored using GNSS RTK.However, reference stations in unstable areas can move, leading to the misleading results.This study proposes a method that uses the PPP technique and CUSUM control chart method to detect the instability of the reference station, and then compensate for the displacements of the monitoring stations.The results show that: The proposed approach offers accurate and trustworthy monitoring results, allowing for precise forecasts and early alerts in the deformation region.The stability analysis of the reference station relies on obtaining accurate reference station coordinates.Due to the lower positioning accuracy and longer convergence time of the PPP technique, the application of the proposed technique to landslide monitoring requires more extensive validation.

( 1 )Fig. 1 Fig. 2
Fig. 1 Flow chart of stability analysis of reference station and compensation for monitoring stations

Fig. 5
Fig. 4 PPP displacement time series (N, E, U) of the reference station TQ01.The data from the first 25 days were selected to evaluate the results of PPP

Fig. 7
Fig. 7 The displacement time series in different directions (N, E, U) of four monitoring stations where TQ02 is used as reference

Table 1
The parameter setting of GNSS software for the Tengqing landslide

Table 4
The cumulative displacement and average rate (N, E, U) of five stations after compensation