A Multi-GNSS Differential Phase Kinematic Post-Processing Method

We propose a multiple global navigation satellite system (multi-GNSS) differential phase kinematic post-processing method, expand the current Track ability, and finely tune processing parameters to achieve the best results for research purposes. The double-difference (DD) phase formulas of GLONASS are especially formulated, and the method of arc ambiguity resolution (AR) in post-processing is developed. To verify the feasibility of this AR method, a group of static baselines with ranges from 8 m to 100 km and two kinematic tests were used. The results imply that 100% of ambiguities in short baselines and over 90% in long baselines can be fixed with the proposed ambiguity resolution method. Better than a 10-mm positioning precision was achieved for all the horizonal components of those selected baselines and the vertical components of the short baselines, and the vertical precision for long baselines is around 20 to 40 mm. In the posterior residual analysis, the means of the residual root-mean-squares (RMSs) of different systems are better than 10 mm for short baselines and at the range of 10–20 mm for baselines longer than 80 km. Mostly, the residuals satisfy the standard normal distribution. It proves that the new method could be applied in bridge displacement and vibration monitoring and for UAV photogrammetry.


Introduction
The global navigation satellite system (GNSS) in recent years has been widely used in many geodetic and geomatic engineering applications. With regard to the kinematic or dynamic applications, such as unmanned aerial vehicle (UAV) measurements [1,2], vehicle-mounted GNSS [3], and structural health monitoring (SHM) [4][5][6], the real-time kinematic (RTK) GNSS method is a primary choice for providing high-precision positioning resolution and, more importantly, in a real-time mode. Even though the double-difference (DD) method can effectively eliminate both of the satellite and receiver clock errors and mitigate the ionospheric and tropospheric delays in the data processing to guarantee the success rate of on-the-fly (OTF) ambiguity fixing, the RTK method cannot always achieve a reliable solution in a harsh environment or when a baseline is longer than 10 km, or even fix wrong ambiguities [7][8][9][10].
GNSS DD kinematic post-positioning is a method to process GNSS data and obtain the single-epoch positioning solutions after data collection [11]. Compared with the data processing in real-time mode, post-positioning could gather more a priori information of the measurements for the benefits of the right

Methods
In GNSS post-processing, the pseudorange measurements are often used to obtain the initial positions and help cycle slip detection and ambiguity resolution. The carrier phase observations are generally used to estimate positions and other parameters. The undifferenced carrier phase observations ϕ can be expressed as [18]: where i and p indicate the satellite and receiver; the subscript m refers to a given frequency; λ m is the wavelength of the frequency m; ρ is the geometric distance; c is the speed of light; and dt p and dt i are the clock errors of the receiver and satellite, respectively. T and I are the slant tropospheric and ionospheric delays, respectively. B i m,p and B i m are the receiver-dependent and satellite-dependent uncalibrated phase delays (UPDs) at the frequency m, respectively. N is the integer ambiguity. ε is the measurement noise of carrier phase observations.
Due to the hardware delays of satellites and receivers, orbit and clock errors, and atmospheric delays, etc., the double-difference method is often applied to eliminate these errors between common-viewing satellites and stations.

GPS/Galileo/BDS DD Observation Equations
We assume that two stations p and q have common-viewing GPS/Galileo/BDS satellites i and j. The DD carrier phase observation equations in the intra-system model for short baselines, for example, less than 10 km, can be written as [18]: where ∇∆ is a DD operator. The residual tropospheric and ionospheric delays are negligible. Thus, only the coordinate parameters and integer ambiguities need to be estimated. In terms of long baselines, which are longer than 10 km, the residual ionospheric and tropospheric errors cannot be neglected. Ionosphere-free (IF) combination observables are used to eliminate the first-order ionospheric delays, and the residual tropospheric errors need to be estimated. The DD IF observations of the carrier phase ∇∆ϕ IF are [18] where λ IF is the wavelength of the IF carrier phase combination. ∇∆N IF denotes the DD IF ambiguity, which can be reformed by an integer wide-lane ambiguity and a narrow-lane ambiguity: where ∇∆N wl = ∇∆N 1 − ∇∆N 2 .

GLONASS DD Observation Equations
GLONASS uses the frequency-division multiple access (FDMA) technique to broadcast navigation signals. The signals in the L1 and L2 bands, which we will call G1 and G2 for the GLONASS dual-frequency carrier phase band in the following statements, are transmitted on 14 adjacent frequencies and the antipodal satellites share the same frequency [18,19]. For the n-th GLONASS satellites, the G1 and G2 carrier frequencies are defined as: where κ k ∈ −7, +6 , and the two-channel frequency separations are ∆ f 1 = 9/16 and ∆ f 2 = 7/16. The basic frequencies of G1 and G2 are f 0 1 = 1602MHz and f 0 2 = 1246MHz, respectively. We should note that ∆ f 1 / f 0 1 = ∆ f 2 / f 0 2 = 1/2848. As shown above, the GLONASS carrier phase observables of different satellites hold different frequencies; therefore, they cannot form DD observations directly as code division multiple access (CDMA) systems. However, we can transform the dual-frequency observations of the carrier phase at original frequencies to common frequencies [20], such as to f 0 1 and f 0 2 as: Thus, the DD model of the GLONASS carrier phase for stations p and q and satellites k and l can be written as: ∇∆ϕ kl m,pq = 2848+κ k 2848 ∆ϕ 0,k m,pq − 2848+κ l 2848 ∆ϕ 0,l m,pq = ∇∆ϕ 0,kl m,pq + κ k 2848 ∆ϕ 0,k m,pq − κ l 2848 ∆ϕ 0,l m,pq Remote Sens. 2020, 12, 2727 4 of 20 According to Sleewagegen et al. [21] and Geng et al. [22], the inter-frequency phase biases are small enough to be neglected for all practical purpose. Therefore, the DD form of GLONASS carrier phase observations in short and long baselines can be expressed as: where ∇∆ N 0,kl m,pq and ∇∆ N 0,kl IF,pq are: The IF combination of ∇∆N 0,kl IF,pq and ∆N 0,l IF,pq has the same formula with (4). As κ l −κ k 2848 is fairly small (from −0.0046 to +0.0046), ∆N 0,l m,pq and ∆N 0,l IF,pq could be estimated in advance and do not have to be fixed. In this case, after ∆N 0,l m,pq and ∆N 0,l IF,pq are applicable, the GLONASS DD ambiguities of (10) and (11) can be estimated and fixed like GPS/Galileo/BDS satellites.

Data Processing Flowchart
Considering the complexity and imperfection of the inter-system difference method, we will only form the GNSS DD observations within intra-systems. The data processing flowchart is shown in Figure 1.
Step 1: Read receiver-independent exchange format (RINEX) observations and broadcast ephemeris, transform the original GLONASS carrier phase observations into the basic frequency as shown in (7), and solve the approximate coordinates of a rover station with DD pseudorange IF combinations (PCs).
Step 2: Detect cycle slips and define ambiguity arcs.
Step 3: Attain initial DD ambiguity for every ambiguity arc. If GLONASS data are available, the SD ambiguities can also be estimated.
Step 4: Form DD observations epoch by epoch, solve the resolution with Kalman filtering (KF) to obtain float solutions, and search and fix ambiguities for ambiguity arcs.
Step 5: Form DD observations, estimate parameters with the fixed ambiguities and use forward and backward KF to smooth the fixed solutions; output the final solutions. The detailed algorithms will be introduced in the following statements.

Data Pre-Processing and Ambiguity Arc Definition
The data pre-processing involves cycle slip detection and short segment cancellation. In the cycle slip detection procedure, the first-order difference of Melbourne-Wübbena (MW) and geometry-free (GF) combination time series will be employed. The MW and GF combinations are formed as: If the difference between adjacent epochs is larger than the threshold (0.25 cycles for GF and 0.5 cycles for MW combinations), an ambiguity flag will be added, and this means an ambiguity arc is ended and another one is started. The new ambiguity arc will be ended at the end of the satellite tracking or another cycle slip epoch. This procedure will be done in the single-differenced observations Remote Sens. 2020, 12, 2727 5 of 20 between stations. Finally, the number of epochs of each ambiguity arc will be counted and the short ones, for example, less than 60 epochs, will be deleted.

Data Pre-Processing and Ambiguity Arc Definition
The data pre-processing involves cycle slip detection and short segment cancellation. In the cycle slip detection procedure, the first-order difference of Melbourne-Wübbena (MW) and geometry-free (GF) combination time series will be employed. The MW and GF combinations are formed as: If the difference between adjacent epochs is larger than the threshold (0.25 cycles for GF and 0.5 cycles for MW combinations), an ambiguity flag will be added, and this means an ambiguity arc is ended and another one is started. The new ambiguity arc will be ended at the end of the satellite tracking or another cycle slip epoch. This procedure will be done in the single-differenced observations between stations. Finally, the number of epochs of each ambiguity arc will be counted and the short ones, for example, less than 60 epochs, will be deleted.

Initial Ambiguity Resolution Using MW+GF Combinations
Before the parameter estimation, we will use MW and GF combinations to calculate the initial ambiguities of the ambiguity arc [23]. The relationship between DD dual-frequency carrier phase ambiguities and DD MW and GF combinations is shown as: Step 1 Step 2 Step 3 Step 4 Step 5 Yes Figure 1. Data processing flowchart.

Initial Ambiguity Resolution Using MW+GF Combinations
Before the parameter estimation, we will use MW and GF combinations to calculate the initial ambiguities of the ambiguity arc [23]. The relationship between DD dual-frequency carrier phase ambiguities and DD MW and GF combinations is shown as: where g = f 2/f 1 . Thus, the DD dual-frequency ambiguity can be calculated by: In the data processing, the DD MW and GF combinations can be calculated only when the overlapped epochs of the common-viewing station and satellite pairs are over 60, and the estimated ambiguities are flagged at the corresponding non-pivot satellite arc of the rover station. We take a weighted average of the MW and GF combinations of each of the overlapped epochs. According to (14), the dual-frequency float ambiguities can be obtained, and we can get the integer initial ambiguities by rounding them directly as: where ∇∆ N 1 and ∇∆ N 2 are the initial integer ambiguities. Then, the residuals of MW and GF combinations, r MW and r GF , after moving the initial integer ambiguities from (13) will be calculated as: Finally, the RMS statistics of the MW and GF combinations of the DD ambiguity arc, σ MW and σ GF , will also be calculated and used in the ambiguity resolution procedure.
It should be noted here that the SD ambiguities of GLONASS arcs need to be estimated in this way as well to make the DD GLONASS ambiguities, which are the integer and estimable in (10) and (11).

Ambiguity Resolution
Before the ambiguity resolution, we firstly estimate the float solutions epoch by epoch using KF. The observation equation of the parameter estimation model can be given as: where the state vector x a contains coordinate parameters of rover stations; x b is the vector for ambiguity parameters and vector x u includes the other parameters. v denotes the noise vector. l is the OMC (observed minus computed observations). P is the weight matrix of DD observations. A, B, and C denote the corresponding coefficient matrices of the parameters, respectively. Thus, the estimated states and their coefficient matrices can be written into vectors x and H as: The parameter states and their covariance matrix will be: where x + e is the a posteriori state vector at epoch e, x − e+1 is the a priori state vector at epoch e+1, Φ e+1,e is the state transition matrix from epoch e to e+1, P + e is the a posteriori covariance matrix of x + e , P − e+1 is the a priori covariance matrix of x − e+1 , and Q e is the process noise covariance matrix. The measurement update step is given as: where K e+1 is the Kalman gain matrix, x + e+1 is the a posteriori parameter state vector at epoch e+1, l e+1 is the measurement residual vector, H e+1 is the measurement design matrix, R e+1 is the measurement noise covariance matrix, and I e+1 is an identity matrix. After the parameter estimation in every epoch, an a posteriori residual editing procedure will be carried out to remove the observations whose residuals are 3 times larger than RMS of DD residuals.
In short baselines, the original L1 and L2 observations will be applied and the estimated parameters are the coordinates of the rover station and dual-frequency ambiguities. Meanwhile, the coordinates of the rover station, IF ambiguities, and residual tropospheric delays will be estimated with IF combinations in long baselines.
Assuming that the estimated float L1 and L2 ambiguities of an ambiguity arc are ∇∆N f ,L1 and ∇∆N f ,L2 with the standard deviations σ L1 and σ L2 in short baselines, we firstly sort all the ambiguities in the order of their standard deviations, and then the ambiguity resolution will be started from the Remote Sens. 2020, 12, 2727 7 of 20 smallest ones. If the float ambiguities and their standard deviations satisfy (21), we can round the float ambiguities as to the fixed ones: where: where relrank defaultly with 25 and σ t with 0.25 cycles are the thresholds.
In long baselines, the estimated float IF ambiguities and standard deviations are ∇∆N f ,IF and σ IF . Due to residual orbit errors, tropospheric delays, and the large noise of IF combinations, we need to search the optimal ambiguities for long baselines. The search space for L1 is from χ min to χ max , and L2 ambiguity is centered on the L1 candidate and takes η as steps on both sides of it, as shown in (23). The ways to estimate χ min , χ max , and η are shown in (24): The candidates will be put into (25) to get the residuals of the IF, MW, and GF combinations. Then, the ambiguity decision process will be based on the residuals to select the best and next best fit values shown in (26) and (27), and is followed by a ratio test using (28). The coefficients a and b in (27) are used to control the effects of MW and GF combinations in the ambiguity resolution. We can adjust a and b to fix more ambiguities and avoid the cases of fixing to wrong ambiguities. Normally, b will be set as 0.1 to reduce the weight of GF combinations in long baselines: If the candidates satisfy all the following conditions, the best candidates will be the fixed ambiguities we searched for. These conditions are: σ IF and σ MW lower than 0.25 and 0.5 cycles, respectively; the best_fit value lower than the threshold with 25.0; the ratio larger than the threshold ξ: Remote Sens. 2020, 12,2727 8 of 20 where iter is the iteration. This means that ξ is gradually decreasing to fix more ambiguities with iterations; however, the correctness of ambiguity fixing will also be checked and controlled by the other conditions. After searching and fixing, the fixed ambiguity(-ies) will be treated as a virtual observation with the weight(s) of 10 9 to update the KF system following (19) to (20), and then, the updated float-valued ambiguity in the following list with the lowest standard deviation will be searched and fixed as the method introduced until the end of the ambiguity list. Then, the KF estimation will be carried out with the fixed ambiguities to estimate the unfixed ones and to search and fix the float-valued ambiguities again. This procedure will be iterated until all the ambiguities are fixed or no more ambiguities can be fixed, which will be kept as a float in the following resolution. In the last step, the data will be resolved with the fixed ambiguities using forward and backward KF to smooth the solutions. Finally, the results will be outputted.

Experimental Analysis
In this section, a group of experiments for collecting static and kinematic data are carried out to evaluate the feasibility of the method proposed above.

Data Collection and Processing Strategies
In the static data collections, we selected baselines of different lengths that GPS, GLONASS, Galileo, and BDS satellites, which can be observed from the Multi-GNSS Experiment (MGEX) network, and the baseline solutions were resolved using GAMIT software in advance to the reference value for comparison purposes. Since not all the receivers can track the BDS-3 satellites, the BDS data applied here were all from the BDS-2 system. The data spanned one week from 24th March 2019 (day of year (DOY) 83) to 30th March 2019 (DOY 89). The kinematic experimental data came from a bridge monitoring platform and a UAV flying trial. The detailed information of the experimental baselines is shown in Table 1. Two observation resolution methods can be selected according to the baseline lengths. For the short baselines, we directly used L1+L2 observations to solve the solutions. Otherwise, IF combinations were used for long baselines. The detailed data processing strategies are listed in Table 2.

Ambiguity Resolution
In this section, we will show the performance of ambiguity resolution. Firstly, the statistics of one day of ambiguity resolution will be given in detail for the data processing strategies with L1+L2 and IF observations. Then, the ambiguity-fixing success rate will be analyzed for all daily sessions. Signals and tracking modes processed The tracking approaches for the bands are sorted in the ascending order of selecting priority, and each tracking mode is represented by one letter [24]: GPS L1/L2: C S L X W GLONASS G1/G2: C P Galileo E1: C X Galileo E5a: I Q X BDS B1/B2: I Q X

Tropospheric parameters
The initial value is zero. The a priori standard deviation is 0.1 m and the variance between epochs is simulated by random walk process, with density 9 × 10 −8 m 2 /s. PCO and PCV GNSS satellites and receiver antennas follow igs14_2029.atx Weighting scheme Elevation dependent model with σ 2 = a 2 + b 2 / sin 2 θ, where θ presents the elevation of satellite.
Ephemeris BRDM combined broadcast ephemeris. Figures 2 and 3 give the detailed statistics of arc ambiguity resolution with L1+L2 and IF observations, respectively, for a whole-day session and the ambiguity arc numbers are listed in the sequence of being fixed. In Figure 2, the standard deviation of the float-valued ambiguities of L1 and L2 are all less than 0.05 cycles, which are significantly lower than the threshold with 0.25 cycles, and the fractional parts of the estimates are confined to within ±0.2 cycles. However, the last few ambiguities show obviously larger values than the previous ones. From the ratio test in the bottom of Figure 2, most of the ratios are much higher than the threshold. The ones with a ratio number of 100 in the figure actually exceed 100, and only 6 out of 144 ambiguities are lower than the threshold. After checking the arc lengths, we know that all of the 6 ambiguity arcs hold only 60-70 epochs, and the elevations are lower than 20 • . Figure 3 shows the conditions that should be satisfied to fix the IF ambiguities in long baselines. The standard deviations of the estimated IF ambiguities and MW combinations are all lower than their corresponding thresholds in the top panel of Figure 3. However, in the middle panel of Figure 3, the fit values of GF combination are relatively larger than IF and MW fit values, especially the last few ambiguities. This is because the GF combinations are significantly affected by the residual ionospheric delays. Hence, we set b in (27) to 0.1 to reduce the weight of the GF combinations in ambiguity resolution. In this case, the fit values of the best fit are all lower than the threshold 25. However, the ratios of 6 ambiguities are all lower than 20 in the iteration, which means they cannot pass the ratio test. By comparing with the fit values, we found that the fittings of these IF ambiguity estimates are also slightly higher than others. This might be due to a poor data quality.

Ambiguity-Fixing Rate
From detailed ambiguity resolution analyses, we list the ambiguity-fixing rate of the daily session for all experiments in Table 3. It shows that over 90% of the ambiguities are fixed in the static data experiments, and it is 100% for short baselines. As for PERT-CUT0, the L1+L2 processing strategy is better than the IF combinations in ambiguity resolution. Probably due to the un-eliminable residual ionospheric and tropospheric delays, the ambiguity resolution rate is lower in the long baseline experiments than that of the short ones. However, it could still achieve over 90% in the experiments.

Static Experimental Analysis
In geophysical research, the kinematic solutions of GNSS permanent stations are applied in general to analyze the daily and sub-daily signals, and further the way they propagate to the annual and semiannual period [28,29]. Besides, the method could also be employed in the monitoring of quasi-static movement of structures, such as dams [30], bridge piers, and geological hazard, etc., to see the spurious signals in the short time scale. Therefore, the positioning precision metrics for different baseline lengths and the corresponding posteriori residual characteristics should be clearly known for an applicable data processing method, software, and engineering application, though we cannot promise the positioning precision in all cases.

Positioning Results
The results of the baseline resolutions after removing the averages are shown in Figure 4. Due to space limitations, we only show the results for DOY 83 and 84 in 2019. In the baseline time series, the positioning precision is getting worse in all three directions with the increase of the baseline length. However, the fluctuation is always within ±2 cm in the horizontal component. The vertical time series change from only 1 cm for CUTA-CUT0 to almost 1 dm for NNOR-PERT and GAMG-DAEJ. Table 4 shows the RMS statistics of the baseline resolutions compared with the GAMIT solutions. It shows that the precisions of the N and E directions are always better than 10 mm for all baselines.   Table 4 shows the RMS statistics of the baseline resolutions compared with the GAMIT solutions. It shows that the precisions of the N and E directions are always better than 10 mm for all baselines.     6 show the residual RMS statistics for every satellite after solving with the L1+L2 observations and IF combinations on DOY 88 for PERT-CUT0. It presents that, mostly, the RMS of satellite residuals is lower than 20 mm. However, they are a little larger than other systems for GLONASS satellites. Table 5 demonstrates the averages of the residual RMS for different satellite systems. Apparently, regardless of the positioning with L1+L2 observations and IF combinations, the means of the residual RMS of different systems are all within 20 mm, better than 10 mm for short baselines and 10-20 mm for long baselines. It should be noted that BDS seems to always have the lowest RMS metrics. From the observation numbers, we know that BDS may take the benefits of more observations of the Geostationary Earth Orbits (GEO) and the Inclined GeoSynchronous Orbit (IGSO) satellite. Figures 5 and 6 show the residual RMS statistics for every satellite after solving with the L1+L2 observations and IF combinations on DOY 88 for PERT-CUT0. It presents that, mostly, the RMS of satellite residuals is lower than 20 mm. However, they are a little larger than other systems for GLONASS satellites. Table 5 demonstrates the averages of the residual RMS for different satellite systems. Apparently, regardless of the positioning with L1+L2 observations and IF combinations, the means of the residual RMS of different systems are all within 20 mm, better than 10 mm for short baselines and 10-20 mm for long baselines. It should be noted that BDS seems to always have the lowest RMS metrics. From the observation numbers, we know that BDS may take the benefits of more observations of the Geostationary Earth Orbits (GEO) and the Inclined GeoSynchronous Orbit (IGSO) satellite.   Figures 5 and 6 show the residual RMS statistics for every satellite after solving with the L1+L2 observations and IF combinations on DOY 88 for PERT-CUT0. It presents that, mostly, the RMS of satellite residuals is lower than 20 mm. However, they are a little larger than other systems for GLONASS satellites. Table 5 demonstrates the averages of the residual RMS for different satellite systems. Apparently, regardless of the positioning with L1+L2 observations and IF combinations, the means of the residual RMS of different systems are all within 20 mm, better than 10 mm for short baselines and 10-20 mm for long baselines. It should be noted that BDS seems to always have the lowest RMS metrics. From the observation numbers, we know that BDS may take the benefits of more observations of the Geostationary Earth Orbits (GEO) and the Inclined GeoSynchronous Orbit (IGSO) satellite.    Figure 7 gives the residual distribution histograms of L1+L2 and IF observations for PERT-CUT0. Basically, over 80% of the residuals of the GPS/BDS/Galileo satellites are within ±2 cm, which is about 0.1 cycles in the carrier phase cycle. The GLONASS one is slightly worse, which is over 70%. Generally, the residuals satisfy the standard normal distribution. For long baselines, however, 80% of GPS and BDS residuals are within ±2 cm, and GLONASS and Galileo are 70%.

Residual Analysis
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 21 Figure 6. The posteriori residual RMS and the observation number statistics of GNSS satellites for PERT-CUT0 after solving with IF combinations on DOY 88.  Figure 7 gives the residual distribution histograms of L1+L2 and IF observations for PERT-CUT0. Basically, over 80% of the residuals of the GPS/BDS/Galileo satellites are within ±2 cm, which is about 0.1 cycles in the carrier phase cycle. The GLONASS one is slightly worse, which is over 70%. Generally, the residuals satisfy the standard normal distribution. For long baselines, however, 80% of GPS and BDS residuals are within ±2 cm, and GLONASS and Galileo are 70%.

Kinematic Experimental Analysis
In this section, two kinematic experiments, including a bridge monitoring experiment and a UAV flying trial, are carried out to show the performance of this method in the practical applications.

Kinematic Experimental Analysis
In this section, two kinematic experiments, including a bridge monitoring experiment and a UAV flying trial, are carried out to show the performance of this method in the practical applications.

Bridge Monitoring Experimental Analysis
The experimental data were from the GeoSHM (GNSS and Earth Observation for Structural Health Monitoring of Bridges) project, which is led by UbiPOS UK Ltd. and participated in by the University of Nottingham, and other partners under sponsorship from the European Space Agency (ESA). The detailed information about GeoSHM can be found in Meng et al. [6,7,31].
In the GeoSHM system, SHM1 is the reference station. SHM2 and SHM3 are two monitoring stations on both sides of the middle span, and SHM4 is the one sitting on the top of the south tower of the Forth Road Bridge (FRB). They are all equipped with an LEICA GM30 receiver and an LEIAR 10 antenna. Figure 8 depicts the locations of stations on the FRB and the baseline lengths. The data were collected on 24 April 2019 UTC 1:00-2:00. Baselines SHM2-SHM1 and SHM4-SHM1 were selected to test the proposed method. Because only one or two BDS satellites can be tracked in the experiments, the BDS observations were out of our consideration in the data processing.
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 21 The experimental data were from the GeoSHM (GNSS and Earth Observation for Structural Health Monitoring of Bridges) project, which is led by UbiPOS UK Ltd. and participated in by the University of Nottingham, and other partners under sponsorship from the European Space Agency (ESA). The detailed information about GeoSHM can be found in Meng et al. [6,7,31].
In the GeoSHM system, SHM1 is the reference station. SHM2 and SHM3 are two monitoring stations on both sides of the middle span, and SHM4 is the one sitting on the top of the south tower of the Forth Road Bridge (FRB). They are all equipped with an LEICA GM30 receiver and an LEIAR 10 antenna. Figure 8 depicts the locations of stations on the FRB and the baseline lengths. The data were collected on 24 April 2019 UTC 1:00-2:00. Baselines SHM2-SHM1 and SHM4-SHM1 were selected to test the proposed method. Because only one or two BDS satellites can be tracked in the experiments, the BDS observations were out of our consideration in the data processing.  Figure 9 shows the displacement time series of the monitoring stations. The Forth Road Bridge mostly approaches from the south to the north. Therefore, we will express the baseline resolutions in the bridge local coordinate frame (NEU). N directs the longitude of the bridge, and E directs the lateral direction and U the vertical direction. In the 3-D displacement time series, SHM2 shows the obvious low-and high-frequency movements in the E direction. This is mainly caused by the wind loading. The time series in the N direction is quite stable and the vibration amplitude is larger in the U direction compared with the horizontal directions. However, the time series of SHM4 is stable since it is located on the top of the supporting tower.   Figure 9 shows the displacement time series of the monitoring stations. The Forth Road Bridge mostly approaches from the south to the north. Therefore, we will express the baseline resolutions in the bridge local coordinate frame (NEU). N directs the longitude of the bridge, and E directs the lateral direction and U the vertical direction. In the 3-D displacement time series, SHM2 shows the obvious low-and high-frequency movements in the E direction. This is mainly caused by the wind loading. The time series in the N direction is quite stable and the vibration amplitude is larger in the U direction compared with the horizontal directions. However, the time series of SHM4 is stable since it is located on the top of the supporting tower.  Figure 10 shows the spectrum analysis of the displacement time series. The remarkable peak frequencies are illustrated by the purple vertical dash lines and are listed in Table 6. Compared with the results from Meng et al. [7], only the first modal frequency of SHM2 in the E direction is different, with 0.065 Hz in Meng et al. [7]. The possible reason is that the vehicle burden has been relieved as only the public transportations are being permitted to pass through the FRB from the beginning of September 2017. Thus, the modal frequency could have changed under different mass loadings due to traffics. The other frequencies showing the same values as Meng et al. indicates that the method proposed is appropriate for the bridge displacement and vibration monitoring application [7].  Figure 10 shows the spectrum analysis of the displacement time series. The remarkable peak frequencies are illustrated by the purple vertical dash lines and are listed in Table 6. Compared with the results from Meng et al. [7], only the first modal frequency of SHM2 in the E direction is different, with 0.065 Hz in Meng et al. [7]. The possible reason is that the vehicle burden has been relieved as only the public transportations are being permitted to pass through the FRB from the beginning of September 2017. Thus, the modal frequency could have changed under different mass loadings due to traffics. The other frequencies showing the same values as Meng et al. indicates that the method proposed is appropriate for the bridge displacement and vibration monitoring application [7].  Figure 10 shows the spectrum analysis of the displacement time series. The remarkable peak frequencies are illustrated by the purple vertical dash lines and are listed in Table 6. Compared with the results from Meng et al. [7], only the first modal frequency of SHM2 in the E direction is different, with 0.065 Hz in Meng et al. [7]. The possible reason is that the vehicle burden has been relieved as only the public transportations are being permitted to pass through the FRB from the beginning of September 2017. Thus, the modal frequency could have changed under different mass loadings due to traffics. The other frequencies showing the same values as Meng et al. indicates that the method proposed is appropriate for the bridge displacement and vibration monitoring application [7].    The residual distribution histograms are shown in Figure 11. It shows that the residuals of all GNSS systems are basically within ±2 cm and satisfy the standard normal distribution.
Remote Sens. 2020, 12, x FOR PEER REVIEW 18 of 21 The residual distribution histograms are shown in Figure 11. It shows that the residuals of all GNSS systems are basically within ±2 cm and satisfy the standard normal distribution. In this experiment, 100% of the ambiguities are fixed for SHM2-SHM1 and 96.2% for SHM4-SHM1 from Table 3. Two factors may affect the ambiguity resolution at SHM4-SHM1, including the obstructions beside SHM4 and the height distance between SHM4 and SHM1, which is over 200 m. The large height distance of the two stations could cause residual tropospheric delays, which may not be completely differenced by the DD method.

UAV Flying Trial Analysis
The UAV flying trial was carried out at Wuhan University in China on 29 March 2019 (DOY 88). It was a UAV photogrammetry trial. The reference station RS01 is located on the roof of the School of Remote Sensing and Information Engineering, which is equipped with a Hi-Target iRTK2 receiver and an IRTK2-2 antenna. The sampling rate was set to 5 Hz and the cutoff elevation is 10°. The UAV is equipped with a DJI PHANTOM 4 RTK equipment. The trial lasted 16 minutes and the flight altitude was about 150 m. The weather was sunny and breezy during the experiment. We collected the GPS, BDS, and GLONASS observations and there were 4684 epochs in total. Because of the equipment failure, only three to four GLONASS satellites could be tracked during the experiment. Therefore, only three GLONASS satellites were available in the data processing. Figure 12 shows the UAV trajectory. In the horizontal plane, a clear aerial photogrammetry type stripes can be observed, and the oscillations in the U direction are shown due to the wind effect. When the UAV changes directions, there would be a larger oscillation shown in the U time series with 10-15 cm; otherwise, the oscillations are mostly within 5 cm. In this experiment, 100% of the ambiguities are fixed for SHM2-SHM1 and 96.2% for SHM4-SHM1 from Table 3. Two factors may affect the ambiguity resolution at SHM4-SHM1, including the obstructions beside SHM4 and the height distance between SHM4 and SHM1, which is over 200 m. The large height distance of the two stations could cause residual tropospheric delays, which may not be completely differenced by the DD method.

UAV Flying Trial Analysis
The UAV flying trial was carried out at Wuhan University in China on 29 March 2019 (DOY 88). It was a UAV photogrammetry trial. The reference station RS01 is located on the roof of the School of Remote Sensing and Information Engineering, which is equipped with a Hi-Target iRTK2 receiver and an IRTK2-2 antenna. The sampling rate was set to 5 Hz and the cutoff elevation is 10 • . The UAV is equipped with a DJI PHANTOM 4 RTK equipment. The trial lasted 16 minutes and the flight altitude was about 150 m. The weather was sunny and breezy during the experiment. We collected the GPS, BDS, and GLONASS observations and there were 4684 epochs in total. Because of the equipment failure, only three to four GLONASS satellites could be tracked during the experiment. Therefore, only three GLONASS satellites were available in the data processing. Figure 12 shows the UAV trajectory. In the horizontal plane, a clear aerial photogrammetry type stripes can be observed, and the oscillations in the U direction are shown due to the wind effect.
When the UAV changes directions, there would be a larger oscillation shown in the U time series with 10-15 cm; otherwise, the oscillations are mostly within 5 cm. Similarly, Figure 13 gives the residual distribution histograms. The residual RMS values of the satellites are all within 10 mm and the residuals are almost 100% within ±2 cm. As for the ambiguity resolution in this experiment, only 71.76% of the ambiguities are fixed in Table 3. It is because the observation tracking mode or channel is inconsistent between the base and rover receivers for several GPS satellites. The DD method cannot eliminate the phase hardware delays, which affects the GPS ambiguity resolution. However, all the BDS ambiguities are successfully fixed, which can guarantee the positioning precision. Similarly, Figure 13 gives the residual distribution histograms. The residual RMS values of the satellites are all within 10 mm and the residuals are almost 100% within ±2 cm. Similarly, Figure 13 gives the residual distribution histograms. The residual RMS values of the satellites are all within 10 mm and the residuals are almost 100% within ±2 cm. As for the ambiguity resolution in this experiment, only 71.76% of the ambiguities are fixed in Table 3. It is because the observation tracking mode or channel is inconsistent between the base and rover receivers for several GPS satellites. The DD method cannot eliminate the phase hardware delays, which affects the GPS ambiguity resolution. However, all the BDS ambiguities are successfully fixed, which can guarantee the positioning precision. As for the ambiguity resolution in this experiment, only 71.76% of the ambiguities are fixed in Table 3. It is because the observation tracking mode or channel is inconsistent between the base and rover receivers for several GPS satellites. The DD method cannot eliminate the phase hardware delays, which affects the GPS ambiguity resolution. However, all the BDS ambiguities are successfully fixed, which can guarantee the positioning precision.

Conclusions
RTK is currently the primary choice in high-precision positioning services. However, it could not always achieve a reliable solution on a harsh environment or even fix wrong ambiguities. A high-precision GNSS data post-processing method is needed to be the back-up solution or augment method to support the numerous geophysical research and geodetic engineering applications. Therefore, we proposed a new multi-GNSS differential phase kinematic post-processing method. We formulated the DD GLONASS phase formulas and demonstrated the arc ambiguity resolution method. Based on the static and kinematic baseline resolutions, the results imply that 100% of the ambiguities in short baselines and over 90% in long baselines can be fixed with the ambiguity resolution method. The positioning precision of better than 10 mm can always be achieved in the horizontal component for the selected baselines, and in the vertical component, with better than 10 mm in short baselines and 30-40 mm in long baselines. In the posterior residual analysis, it was found that the RMS of the residuals for different systems were all better than 10 mm. Mostly, the residuals satisfy the standard normal distribution. Besides, we demonstrated that this method could be used for bridge displacement and vibration monitoring and UAV photogrammetry.