Application of Multi-System Combination Precise Point Positioning in Landslide Monitoring

To verify the positioning performance and reliability of multi-system combination Precise Point Positioning in landslide monitoring, we carried out a multi-system combination Precise Point Positioning calculation experiment on the monitoring data of a single landslide disaster area in Fujian Province. The coordinates of the monitoring points obtained by a continuously operating reference station and the monitoring station for static relative positioning were used as reference values. The GPS system was used as the standard system and the combined PPP solution mode of G/R/C, G/R/E and G/R/E/C was used to obtain the surface displacement of the landslide area. The research showed that multi-system combination PPP converges to the centimeter level in about 30 min. The average value of internal accordant precision was more than 1 mm after convergence, and that of the external accordant precision was more than 5 cm, which meets the centimeter-level accuracy requirements in rapid landslide deformation monitoring.


Introduction
As a sudden geological disaster, landslides are highly destructive and severely harmful to people's lives, property and safety. As the source of disasters is usually located in the middle and upper part of mountains and the vegetation coverage is high in most areas, it is often difficult for geological survey personnel to reach them, and thus these areas are characterized by strong concealment and difficult investigation [1,2]. The questions of how to strengthen the monitoring of the hidden danger of geological disaster, make timely discoveries and effectively identify the potential landslide disaster have become the focus and difficulty of geological disaster prevention and control.
Traditional landslide deformation monitoring uses a total station and level to monitor the displacement changes in the slope development process, which has some shortcomings, such as long monitoring periods, the significant influence of terrain, difficult measurement, high labor intensity and poor real-time performance [3][4][5]. Since the 1990s, differential GNSS positioning technology has been widely used in geological hazard monitoring [6,7]. By using two GNSS receivers to form double difference observation values, common errors such as the receiver clock error and satellite clock error can be eliminated, and strong correlation errors such as tropospheric delay and ionospheric delay can be lessened. The positioning accuracy can reach centimeter level or even millimeter level [8,9]. However, differential GNSS positioning technology requires at least one receiver to be placed on a known station for observation. When the number of visible satellites is limited, the distance between the known point and the monitoring point is too far, the environment is complex, and the height difference is too large, then the positioning accuracy is not high [10]. In 1997, American Jet Power Laboratory (JPL) [11] presented a high-precision, non-differential GPS single-point absolute positioning technology called Precise Point Positioning (PPP). It is made of a single GNSS receiver, using an IGS (International GNSS Services) center to provide the precision of satellite orbit and clock products, for a variety of model to intensive model error corrections. It uses the acquisition of pseudo-range and phase observation values to obtain high precision for coordinates of the designated position according to the earth's international frame of reference. Without the support of the reference station, the accuracy of PPP can reach a decimeter or centimeter scale. Under the condition of long-time static observation, it can achieve the same accuracy as traditional differential GPS positioning [12]. Because of the rapid development of the multi-system GNSS and the growing maturity of precision satellite orbit and precision clock calibration products provided by IGS, the application of multi-system GNSS PPP technology is the current research hotspot. Domestic and foreign scholars have studied GPS/BDS [13][14][15][16], GPS/GLONASS [17][18][19], BDS/GPS/GLONASS [20,21], GPS/BDS/Galileo/GLONASS [22][23][24][25] and other combination modes. A large number of studies have shown that the PPP results of different combinations are different in terms of convergence time and positioning accuracy. In the case of fewer visible satellites and a complex surrounding environment, multi-system combination PPP has obvious advantages. In the application of PPP technology in landslide monitoring, Wang [26] used GIPSY's PPP module to calculate the GPS data of peristaltic landslides for two years in a row. The study showed that the resolution accuracy of horizontal and elevation directions with an observation duration of 24 h was 2-3 mm and 8 mm, respectively. When there was no extreme weather in the observation period, the accuracy of the horizontal direction was 5 mm for the observation period of 4 h, and the accuracy of the elevation direction was 10 mm for the observation period of 8 h. Capillia et al. [27] used a special instrument to carry out a landslide simulation experiment and used PPP technology to solve the problem. The study showed that the average deviation of real-time PPP in the north, east and elevation directions was 2 cm. Zhang et al. [28] calculated 6-12 hours of monitoring data of a landslide area in southern Gansu through GAMIT and PANDA PPP modules. The experimental results showed that using the GAMIT PPP module in calculating the coordinates of the accuracy and precision of relative positioning algorithms to obtain the coordinates of the amount of change can be reflected to five decimal places, and GAMIT and PANDA PPP modules of coordinate difference variation can be reflected to more than six decimal places, which satisfies the precision requirement of rapid landslide monitoring. Huang et al. [29] discussed the calculation accuracy of GNSS PPP and RTK in a large-scale loess landslide, and the study showed that the BDS/GPS combination significantly improved the accuracy compared with the single system. The positioning accuracy of RTK was better than 3 mm and could be used for real-time monitoring of the landslide, while the positioning accuracy of PPP was better than 6 cm and could be used for the maintenance of deformation data. Peng et al. [30] took the data of eight monitoring stations in the Xishan Village landslide area of Sichuan Province for 16 months as the data source and discussed the positioning accuracy and reliability of GPS PPP, BDS PPP and GPS/BDS combined PPP. The results showed that the mean square error of GPS PPP, GPS/BDS PPP, GPS/BDS combined PPP was about 1 cm, 2 cm, and 5 cm, respectively, which meets the accuracy requirements of landslide monitoring in mountainous areas.
At present, the application of PPP in landslide monitoring is mostly single-system PPP or dual-system combined PPP [31,32], and scholars have barely discussed the accuracy, reliability and calculation rate of multi-system combined PPP in landslide monitoring in a continuous period. Their conclusions show that the combined system has greater improvement in positioning stability, continuity and positioning accuracy than the single positioning system, but the specific values are not shown. Therefore, based on the static experiment and dynamic simulation experiment research, this paper experimented with monitoring data of one single landslide disaster area in Fujian province. The coordinates of the monitoring points obtained by static relative positioning between the continuously running reference station and the monitoring station are used as the reference value, and the GPS system is used as the reference system, showing the multi-system combination of G/R/C, G/R/E and G/R/E/C to have better positioning performance, which verifies the positioning performance and reliability of multi-system combination PPP in landslide monitoring.

PPP Model and Evaluation Index
The traditional PPP function model is a combination of eliminating the ionosphere between the dual-frequency pseudo-distance and the phase observation. When a GNSS precision product is used, the observation equation is [33,34]: where: S = satellite system; IF = ionospheric cancellation combination; P = observation value of pseudo-distance; L = phase observation; ρ = geometric distance; c = speed of light; T S = tropospheric delay; dt S = receiver clock difference; λ S = wavelength; N S = ambiguity; ε P S IF = pseudo-range observation noise; ε L S IF = phase observation noise. The error caused by the earth's tides, phase entanglement and earth's rotation have been corrected according to the corresponding model and therefore are not listed in the formula.
Inter-system Bias (ISB) parameter estimation is required when multi-system GNSS combined PPP is carried out. The observation equation of dual-system combined PPP is: where: B = reference satellite system; S i = rest of the constellations except the base system, i = 1,2,3; IBS B,S i = system deviation between the base system and the rest of the system.
As the tropospheric error, ionospheric error and multipath error are strongly correlated with the satellite altitude angle, in order to reduce the impact of errors on PPP positioning performance, the altitude angle stochastic model is adopted in this paper. As the noise and multi-path errors of the satellite observations are concentrated in the low altitude satellite, in order to reduce the weight of the low altitude satellite observations, the piecewise weighting strategy is adopted in this paper. The calculation formula of the pseudorange observation and phase observation is as follows [35]: where: σ 2 P,0 = prior variances of pseudo-range; σ 2 L,0 = phase observations; α = the height of the angle threshold (30 • ).
The signal of the GLONASS system is transmitted by frequency-division multiple access, which leads to different hardware delays in its receiving channels. If the hardware delay is estimated as an unknown parameter, the parameter estimation result will be unstable due to too many parameters being estimated. Therefore, the processing strategy of giving a small weight to the observed value of GLONASS is usually adopted [18]. Therefore, the empirical model is used to set the ratio of pseudo-distance standard deviation and phase standard deviation of the reference system and auxiliary system: σ B,P : σ S 1 ,P : σ S 2 ,P : σ S 3 ,P = 1 : 2 : 2 : 2 σ B,L : σ S 1 ,L : σ S 2 ,L : σ S 3 ,L = 1 : 2 : 2 : 2 (4) where: σ B,P = standard deviation of pseudo-range observations of reference system (3 m); σ S i ,P = standard deviation of pseudo-range observations of auxiliary system (6 m); σ B,L = standard deviation of phase observation of reference system (0.003 m); σ S i ,L = standard deviation of phase observation of auxiliary system (0.006 m).

Convergence Time
Convergence time refers to the number of periods required from the first period to the time when the calculation accuracy meets the limit requirements, which includes the threshold value of point-position accuracy and the number of continuous periods. If an epoch of one-dimensional (north (N), east (E), elevation (U)), two-dimensional (2D) or three-dimensional (3D) is within the scope of the deviation between the set point precision where: CT = convergence time (min); SI = sampling interval; L = the number of epochs required from the beginning of the first epoch until the calculation accuracy meets the limitation requirements.

Root Mean Square Error
The Root Mean Square Error (RMSE) reflects the difference between the estimated value and the reference value, and is expressed as: where: k = initial epoch at the time of convergence time; n = total epoch; P i = parameter to be estimated; P = the reference value.

Overview of the Study Area
The study area is located in the northwest of Quanzhou City, Fujian Province ( Figure 1). It is a medium-low mountain landform of tectonic erosion, with the highest elevation of 957 m, the lowest elevation of 290 m and the height difference of 650 m. The landslide area is located at the foot of the middle and front slope. Topographically, the boundary of 470 m and 350 m constitutes the junction of steep and slow slopes. The elevation of 290 m~350 m is the gentle slope, and the slope is about 5 •~1 0 • . The elevation of 350 m~470 m is the area where villagers gather in Yaoshan. The human engineering reconstruction is intense. The overall slope is about 15 •~2 0 • , and the slope of some artificial slopes behind the house can reach 50 •~7 0 • . The elevation of 470 m~957 m is a steep slope, and the slope is about 35 •~4 0 • . Slope body planting of tea, vegetables and other cash crops, and dense vegetation, has a coverage rate of 70~80%.

Root Mean Square Error
The Root Mean Square Error (RMSE) reflects the difference between the estimated value and the reference value, and is expressed as: where: k = initial epoch at the time of convergence time; n = total epoch; Pi = parameter to be estimated; P = the reference value.

Overview of the Study Area
The study area is located in the northwest of Quanzhou City, Fujian Province ( Figure  1). It is a medium-low mountain landform of tectonic erosion, with the highest elevation of 957 m, the lowest elevation of 290 m and the height difference of 650 m. The landslide area is located at the foot of the middle and front slope. Topographically, the boundary of 470 m and 350 m constitutes the junction of steep and slow slopes. The elevation of 290 m~350 m is the gentle slope, and the slope is about 5°~10°. The elevation of 350 m~470 m is the area where villagers gather in Yaoshan. The human engineering reconstruction is intense. The overall slope is about 15°~20°, and the slope of some artificial slopes behind the house can reach 50°~70°. The elevation of 470 m~957 m is a steep slope, and the slope is about 35°~40°. Slope body planting of tea, vegetables and other cash crops, and dense vegetation, has a coverage rate of 70~80%.

Basic Characteristics of Landslides in the Study Area
There are many transverse and oblique tensile cracks in the slope in the study area and they are mainly scattered highway cracks, up to dozens of which are distributed in the elevation of 360 m~470 m. They are perpendicular to the direction of the highway and intersect locally. The shape is zigzag, the nature is mainly subsidence-type tension, and

Basic Characteristics of Landslides in the Study Area
There are many transverse and oblique tensile cracks in the slope in the study area and they are mainly scattered highway cracks, up to dozens of which are distributed in the elevation of 360 m~470 m. They are perpendicular to the direction of the highway and intersect locally. The shape is zigzag, the nature is mainly subsidence-type tension, and locally it is shearing. The extension length is about 5 m~15 m, the opening is about 5 cm~15 cm, and the visible depth is 10 cm~20 cm. The rock mass near the air side near the local crack also shows the characteristics of down-fault, with a height of about 10 cm~20 cm (Figures 2  and 3). The overall characteristics of the slope body are "heavy rain big slippery, light rain small slippery, no rain not slippery". locally it is shearing. The extension length is about 5 m~15 m, the opening is about 5 cm~15 cm, and the visible depth is 10 cm~20 cm. The rock mass near the air side near the local crack also shows the characteristics of down-fault, with a height of about 10 cm~20 cm (Figures 2 and 3). The overall characteristics of the slope body are "heavy rain big slippery, light rain small slippery, no rain not slippery".

Location of Monitoring Points
According to a large amount of survey data in the early stage, the deformation range of the landslide study area, the position of the main sliding surface and the position of the monitoring point are determined (Figure 4). Within the solid red line is the influence range of the landslide area. The location of GNSS monitoring points ( Figure 5) is near the central axis of the landslide study area, which can effectively monitor the displacement variation of the monitoring points in the study area in real time. The continuously operating monitoring station is equipped with an N72 receiver, A220GR choke antenna, solar panels, lightning rods and other equipment. locally it is shearing. The extension length is about 5 m~15 m, the opening is about 5 cm~15 cm, and the visible depth is 10 cm~20 cm. The rock mass near the air side near the local crack also shows the characteristics of down-fault, with a height of about 10 cm~20 cm (Figures 2 and 3). The overall characteristics of the slope body are "heavy rain big slippery, light rain small slippery, no rain not slippery".

Location of Monitoring Points
According to a large amount of survey data in the early stage, the deformation range of the landslide study area, the position of the main sliding surface and the position of the monitoring point are determined (Figure 4). Within the solid red line is the influence range of the landslide area. The location of GNSS monitoring points ( Figure 5) is near the central axis of the landslide study area, which can effectively monitor the displacement variation of the monitoring points in the study area in real time. The continuously operating monitoring station is equipped with an N72 receiver, A220GR choke antenna, solar panels, lightning rods and other equipment.

Location of Monitoring Points
According to a large amount of survey data in the early stage, the deformation range of the landslide study area, the position of the main sliding surface and the position of the monitoring point are determined (Figure 4). Within the solid red line is the influence range of the landslide area. The location of GNSS monitoring points ( Figure 5) is near the central axis of the landslide study area, which can effectively monitor the displacement variation of the monitoring points in the study area in real time. The continuously operating monitoring station is equipped with an N72 receiver, A220GR choke antenna, solar panels, lightning rods and other equipment.

Application of Multi-system Combination PPP in Landslide Monitoring
Using the landslide monitoring data of 18 June 2019~18 January 2020 in the research area, the monitoring data were processed and analyzed by static relative positioning technology, and the displacement of the monitoring point was obtained. The relative accuracy of the baseline vector is better than 0.1 ppm, the position accuracy of the plane direction is better than 2 mm and the position accuracy of the elevation direction is better than 5 mm. The traditional measurement methods through precision level measurement validate the results of the static relative positioning. The results show that the static relative positioning monitoring results and precision level measurement results are basically identical and have high precision and reliability. In this paper, the coordinates of monitoring points obtained by static relative positioning technology are taken as reference values, and the convergence time, positioning accuracy and stability of multi-system combined PPP technology in landslide monitoring are emphatically discussed. To be concise, GPS, GLONASS, Galileo and BDS are represented by the letters G, R, E and C, respectively.

Analysis of Convergence Time
Combined with the data quality test results of monitoring points and the preliminary experimental conclusions, the GPS system was used as the benchmark to carry out the multi-system combined PPP calculation experiment. It used G/R/C, G/R/E and G/R/E/C as three good combinations and the static relative positioning calculation coordinate value as the reference coordinate to verify the positioning performance and reliability of G/R/C, G/R/E and G/R/E/C combined PPP. As the calculated results of the observed data are basically similar, the position deviation and observation status of the observed data on the

Application of Multi-system Combination PPP in Landslide Monitoring
Using the landslide monitoring data of 18 June 2019~18 January 2020 in the research area, the monitoring data were processed and analyzed by static relative positioning technology, and the displacement of the monitoring point was obtained. The relative accuracy of the baseline vector is better than 0.1 ppm, the position accuracy of the plane direction is better than 2 mm and the position accuracy of the elevation direction is better than 5 mm. The traditional measurement methods through precision level measurement validate the results of the static relative positioning. The results show that the static relative positioning monitoring results and precision level measurement results are basically identical and have high precision and reliability. In this paper, the coordinates of monitoring points obtained by static relative positioning technology are taken as reference values, and the convergence time, positioning accuracy and stability of multi-system combined PPP technology in landslide monitoring are emphatically discussed. To be concise, GPS, GLONASS, Galileo and BDS are represented by the letters G, R, E and C, respectively.

Analysis of Convergence Time
Combined with the data quality test results of monitoring points and the preliminary experimental conclusions, the GPS system was used as the benchmark to carry out the multi-system combined PPP calculation experiment. It used G/R/C, G/R/E and G/R/E/C as three good combinations and the static relative positioning calculation coordinate value as the reference coordinate to verify the positioning performance and reliability of G/R/C, G/R/E and G/R/E/C combined PPP. As the calculated results of the observed data are basically similar, the position deviation and observation status of the observed data on the

Application of Multi-System Combination PPP in Landslide Monitoring
Using the landslide monitoring data of 18 June 2019~18 January 2020 in the research area, the monitoring data were processed and analyzed by static relative positioning technology, and the displacement of the monitoring point was obtained. The relative accuracy of the baseline vector is better than 0.1 ppm, the position accuracy of the plane direction is better than 2 mm and the position accuracy of the elevation direction is better than 5 mm. The traditional measurement methods through precision level measurement validate the results of the static relative positioning. The results show that the static relative positioning monitoring results and precision level measurement results are basically identical and have high precision and reliability. In this paper, the coordinates of monitoring points obtained by static relative positioning technology are taken as reference values, and the convergence time, positioning accuracy and stability of multi-system combined PPP technology in landslide monitoring are emphatically discussed. To be concise, GPS, GLONASS, Galileo and BDS are represented by the letters G, R, E and C, respectively.

Analysis of Convergence Time
Combined with the data quality test results of monitoring points and the preliminary experimental conclusions, the GPS system was used as the benchmark to carry out the multi-system combined PPP calculation experiment. It used G/R/C, G/R/E and G/R/E/C as three good combinations and the static relative positioning calculation coordinate value as the reference coordinate to verify the positioning performance and reliability of G/R/C, G/R/E and G/R/E/C combined PPP. As the calculated results of the observed data are basically similar, the position deviation and observation status of the observed data on the third day in the north (N), east (E) and elevation (U) directions are taken as an example (Figures 6-8).
Appl. Sci. 2021, 11, x FOR PEER REVIEW 8 of 15 third day in the north (N), east (E) and elevation (U) directions are taken as an example (Figures 6-8).       According to the analysis of Figures 6-8, the G/R/C combination observation state is good, the average number of available satellites is 20.6, the average PDOP value is 1.36, and the three directions of NEU converge to within 0.1 m after 12.5, 14 and 30.5 min, respectively. After convergence, the positioning accuracy of N direction is better than 1 cm, E direction is better than 2 cm, and U direction is better than 5 cm. The observation state of the G/R/E combination is slightly worse: the average number of available satellites is 12.4, and the average PDOP value is 1.79. After 12.5, 13.5 and 33 min, the three directions of NEU converge to within 0.1 m. After convergence, the positioning accuracy of N direction is about 1 cm, E direction is better than 2.2 cm, and U direction is better than 5 cm. Because the station is located in the Fujian province, Galileo has fewer visible satellites, G/R/E/C composite state of observation is similar to G/R/C composite, the average number of available satellites is 21.1, and the average PDOP value of 1.35. The three directions of NEU converge to within 0.1 m after 12.5, 14 and 30.5 min, respectively. After convergence, the positioning accuracy of N direction is better than 1 cm, E direction is better than 2 cm, and U direction is better than 5 cm. In summary, the PPP convergence time of multi-system combination is basically the same, and the positioning accuracy of NEU in all three directions can reach centimeter level after about 30 min.
In the combined PPP solution of multi-day GNSS monitoring data, if data stitching is not carried out, the daily observation data will experience about 30 min of initialization time, causing a serious waste of the original observation data. Therefore, it is extremely important to mosaic the multi-day data. In this study, the CHCData Mosaic software provided by Huace Navigation was used to mosaic the original observation data of the monitoring point, and the multi-system combined PPP calculation experiment was carried out, which greatly improved the utilization rate of the original observation. Figure 9 shows the solution results of the observation data stitching on the third and fourth days of the G/R/C, G/R/E and G/R/E/C combinations, respectively. The solid black line at the 2880 epoch is the splice point of observation data in the two phases. According to the analysis of Figure  9, the positioning accuracy of the combination of G/R/C and G/R/E/C in the N direction is According to the analysis of Figures 6-8, the G/R/C combination observation state is good, the average number of available satellites is 20.6, the average PDOP value is 1.36, and the three directions of NEU converge to within 0.1 m after 12.5, 14 and 30.5 min, respectively. After convergence, the positioning accuracy of N direction is better than 1 cm, E direction is better than 2 cm, and U direction is better than 5 cm. The observation state of the G/R/E combination is slightly worse: the average number of available satellites is 12.4, and the average PDOP value is 1.79. After 12.5, 13.5 and 33 min, the three directions of NEU converge to within 0.1 m. After convergence, the positioning accuracy of N direction is about 1 cm, E direction is better than 2.2 cm, and U direction is better than 5 cm. Because the station is located in the Fujian province, Galileo has fewer visible satellites, G/R/E/C composite state of observation is similar to G/R/C composite, the average number of available satellites is 21.1, and the average PDOP value of 1.35. The three directions of NEU converge to within 0.1 m after 12.5, 14 and 30.5 min, respectively. After convergence, the positioning accuracy of N direction is better than 1 cm, E direction is better than 2 cm, and U direction is better than 5 cm. In summary, the PPP convergence time of multisystem combination is basically the same, and the positioning accuracy of NEU in all three directions can reach centimeter level after about 30 min.
In the combined PPP solution of multi-day GNSS monitoring data, if data stitching is not carried out, the daily observation data will experience about 30 min of initialization time, causing a serious waste of the original observation data. Therefore, it is extremely important to mosaic the multi-day data. In this study, the CHCData Mosaic software provided by Huace Navigation was used to mosaic the original observation data of the monitoring point, and the multi-system combined PPP calculation experiment was carried out, which greatly improved the utilization rate of the original observation. Figure 9 shows the solution results of the observation data stitching on the third and fourth days of the G/R/C, G/R/E and G/R/E/C combinations, respectively. The solid black line at the 2880 epoch is the splice point of observation data in the two phases. According to the analysis of Figure 9, the positioning accuracy of the combination of G/R/C and G/R/E/C in the N direction is better than 1 cm, the E direction is better than 2 cm, and the U direction is better than 5 cm after the PPP convergence. After the convergence of the G/R/E combination PPP, the positioning accuracy in the N direction is about 1 cm, the positioning accuracy in the E direction is better than 2.2 cm, and the positioning accuracy in the U direction is better than 5 cm.
better than 1 cm, the E direction is better than 2 cm, and the U direction is better than 5 cm after the PPP convergence. After the convergence of the G/R/E combination PPP, the positioning accuracy in the N direction is about 1 cm, the positioning accuracy in the E direction is better than 2.2 cm, and the positioning accuracy in the U direction is better than 5 cm.

Analysis of Positioning Accuracy
To further verify the reliability of the calculating results, continuous reference stations were set up at the stable points outside the landslide disaster area for static synchronous observation. The calculation coordinates of 24 h static synchronous observation data between the reference point and the monitoring point were taken as the initial coordinates of the monitoring point. The HCMonitor software provided by Huace Navigation was used to calculate the static relative positioning. Figure 10 shows the cumulative changes of the static relative positioning of the monitoring points in the north (N), east (E) and elevation (U) directions. According to the analysis of Figure 10, the cumulative displacement of the monitoring point data in the observation period in the N direction is 9.54 mm, the cumulative displacement in the E direction is 7.38 mm, and in the U direction, it is 19.2 mm. During the observation period, the horizontal displacement of the monitoring point moved along the N and E directions, with the trend of slightly sinking in the elevation direction. The static relative positioning results can identify the millimeter-scale landslide shape variables, which provides a reliable external verification means for PPP results.

Analysis of Positioning Accuracy
To further verify the reliability of the calculating results, continuous reference stations were set up at the stable points outside the landslide disaster area for static synchronous observation. The calculation coordinates of 24 h static synchronous observation data between the reference point and the monitoring point were taken as the initial coordinates of the monitoring point. The HCMonitor software provided by Huace Navigation was used to calculate the static relative positioning. Figure 10 shows the cumulative changes of the static relative positioning of the monitoring points in the north (N), east (E) and elevation (U) directions. According to the analysis of Figure 10, the cumulative displacement of the monitoring point data in the observation period in the N direction is 9.54 mm, the cumulative displacement in the E direction is 7.38 mm, and in the U direction, it is 19.2 mm. During the observation period, the horizontal displacement of the monitoring point moved along the N and E directions, with the trend of slightly sinking in the elevation direction. The static relative positioning results can identify the millimeter-scale landslide shape variables, which provides a reliable external verification means for PPP results.   Table 1 quantifies the difference between the results of the PPP solution and the static relative positioning solution in the three directions of NEU at an interval of one month. To sum up, the accuracy of the calculation results of G/R/C, G/R/E and G/R/E/C combined PPP is basically the same, and the external coincidence accuracy is better than 5 cm for all solutions, so they can be applied to the landslide monitoring of centimeter-level slope movement.   Table 1 quantifies the difference between the results of the PPP solution and the static relative positioning solution in the three directions of NEU at an interval of one month. To sum up, the accuracy of the calculation results of G/R/C, G/R/E and G/R/E/C combined PPP is basically the same, and the external coincidence accuracy is better than 5 cm for all solutions, so they can be applied to the landslide monitoring of centimeter-level slope movement.

Discussion
(1) In the early stage of the experiment, we conducted a single-system, dual-system and system combined static PPP calculating experiment to determine the benchmark system to carry on the static simulation dynamic experiment with a better applied combination in the measured dynamic data. It has important practical significance and application value for the analysis of current positioning performance of the combined system dynamic PPP.
(2) Early experiments on the measured dynamic data decoding showed that in the measured dynamic data quality that the IGS stations provide, the static data quality is poor. This is because the measured dynamic data observation environment is more complex, satellites change frequently, cycles slip and gross error is more often caused. In the future, with the BDS and Galileo precision products, precision ascension and receiver-end and PCV PCO model elaboration, the combined system dynamic PPP positioning accuracy and convergence time will improve further.
(3) For the combined system presented in this paper, in the PPP calculating monitoring data we only used the rapid ephemeris products offered by GFZ. In future research, this can be combined with the CNES, BKG and ESA real-time precision products offered by institutions to determine the accuracy of the evaluation of landslides. Furthermore, it can also be combined with the relevant data such as the comprehensive meteorological, hydrological, geological evaluation of the stability of landslides and deep discussion on the mechanism of the formation of landslides.

Conclusions
In this paper, the coordinates of monitoring points obtained by static relative positioning of the continuously operating reference station and monitoring station were taken as reference values. The GPS system was taken as the reference system, and G/R/E, G/R/C and G/R/E/C combined PPP calculation modes were used to carry out multi-system combined PPP calculation experiments on the data of monitoring points. The results show that: (1) G/R/C, G/R/E, and G/R/E/C combined PPP converged to within 0.1 m after 12.5, 14 and 30.5; 12.5, 13.5 and 33; and 12.5, 14 and 30.5 minutes in NEU directions, respectively. The PPP convergence time of the three combinations is basically the same, and the positioning accuracy of all the three directions of NEU can reach centimeter level after about 30 min.
(2) The accuracy of the PPP solution results of G/R/C, G/R/E and G/R/E/C combinations is basically the same. After the three PPP combinations converge in the NEU directions, the internal coincidence accuracy of the positioning results is better than 1 mm on average, and the external coincidence accuracy is better than 5 cm, which is relatively stable and suitable for landslide monitoring with the slope movement of a centimeter.

Data Availability Statement:
The multi-GNSS precise orbit and clock products issued provided by IGS are available from ftp://igs.ensg.ign.fr/pub/igs/products/ (accessed on 30 June 2021).