A Multi-GNSS Differential Phase Kinematic Post-Processing Method

Abstract

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.

1. 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 ambiguity resolution. At the same time, the quality of pseudorange and carrier-phase observations can be evaluated in advance to get an accurate weighting scheme and remove bad observations. Therefore, the DD kinematic post-positioning method could achieve a higher ambiguity fixing rate and better positioning precision [11].

One of the most famous DD kinematic post-positioning software is Track module in the GAMIT/GLOBK (GNSS At MIT and GLOBal Kalman filter) package. Track supports multi-baseline global positioning system (GPS) data processing in the post-mode. The principle of Track is to resolve and fix the ambiguities first, then estimate the positioning parameters of every epoch using the Kalman filter (KF), and finally attain the trajectory of the moving carrier [12]. The positioning precision of Track is baseline-length dependent, with a few millimeters to the centimeter level [13]. Track has been gradually evolved in geophysical and geodetic research, such as earthquake monitoring by SHM et al. [14,15,16]. However, it can only process GPS data and relevant research and applications have not been explored yet.

Currently, there are four GNSSs available for users, GPS, BeiDou Navigation Satellite System (BDS), ГЛOбaльнaя HAвигaциoннaя Cпyтникoвaя Cиcтeмa (GLONASS), and Galileo Navigation Satellite System (Galileo). Since the last BDS-3 satellite was launched into the designated orbit on 23th June 2020, the BDS-3 system has been completely constructed. Therefore, there are three GNSSs fully available to provide precise positioning, navigation, and timing (PNT) services, and Galileo is to become fully operational at the end of 2020/beginning of 2021. In the future, more than 100 GNSS satellites will be available and 35–40 satellites (depending on the mask angel) are visible for the user in areas with no obstructions. Thus, the integrated multiple GNSS (multi-GNSS) processing will be the trend and will be utilized in positioning, navigation, and timing applications. The multi-GNSS-based precise point positioning (PPP) and RTK algorithms have also been extensively studied in recent years [17].

However, the DD kinematic post-positioning methods with multi-GNSS observations have rarely been investigated. It is important as it can be used as a back-up solution or augment method to support the precise positioning service in cases when the real-time solutions are either not available or not stable. Therefore, we expand the current Track software to support data processing and finely tune the processing parameters to achieve the best results for research purposes. The detailed DD models of the multi-GNSS observations, especially the phase DD model of GLONASS, are formulated and the segmental arc ambiguity resolution method in a post-processing mode is demonstrated. Following several experiments, the ambiguity-fixing performance, positioning precision statistics, and posterior residual analysis is carried out to evaluate the performance of the proposed method. 2. 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]:

_{λm φm,pi}=_ρpi+c(d_tp−d^ti)+_Tm,pi−_Im,pi+_λm(_Nm,pi+_Bm,pi−_Bmi)+_εm,pi,

whereiandpindicate the satellite and receiver; the subscriptmrefers to a given frequency;_λmis the wavelength of the frequencym;ρis the geometric distance;cis the speed of light; andd_tpandd^tiare the clock errors of the receiver and satellite, respectively.TandIare the slant tropospheric and ionospheric delays, respectively._Bm,piand_Bmiare the receiver-dependent and satellite-dependent uncalibrated phase delays (UPDs) at the frequency m, respectively.Nis 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. 2.1. GPS/Galileo/BDS DD Observation Equations

We assume that two stationspandqhave common-viewing GPS/Galileo/BDS satellitesiandj . 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]:

_λm∇Δ_φm,pqij=∇Δ_ρpqij+_λm∇Δ_Nm,pqij+∇Δ_εm,pqij,

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]

_λIF∇Δ_φIF,pqij=∇Δ_ρpqij+∇Δ_Tpqij+_λIF∇Δ_NIF,pqij+∇Δ_εIF,pqij,

where_λIFis the wavelength of the IF carrier phase combination.∇Δ_NIFdenotes the DD IF ambiguity, which can be reformed by an integer wide-lane ambiguity and a narrow-lane ambiguity:

∇Δ_NIF=α·∇Δ_N1−β·∇Δ_N2=γ·∇Δ_N1+β·∇Δ_Nwlα=_f12f12−_f22,β=_{f1 f2f12}−_f22,γ=_f1f1+_f2

where∇Δ_Nwl=∇Δ_N1−∇Δ_N2.

2.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:

_fmk=_fm0+^κkΔ_fm,m=1,2,k=1,⋯,n,

where^κk∈⌈−7,+6⌉, and the two-channel frequency separations areΔ_f1=9/16andΔ_f2=7/16. The basic frequencies of G1 and G2 are_f10=1602MHzand_f20=1246MHz, respectively. We should note thatΔ_f1/_f10=Δ_f2/_f20=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_f10and_f20as:

_φm,pk=_fm0+^κkΔ_{fmfm0φm,p0,k}=2848+^κk2848_φm,p0,k.

Thus, the DD model of the GLONASS carrier phase for stations p and q and satellites k and l can be written as:

∇Δ_φm,pqkl=2848+^κk2848Δ_φm,pq0,k−2848+^κl2848Δ_φm,pq0,l=∇Δ_φm,pq0,kl+(^κk2848Δ_φm,pq0,k−^κl2848Δ_φm,pq0,l)

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:

_λmkΔ_φm,pqk−_λmlΔ_φm,pql=∇Δ_ρpqkl+_λm0∇Δ_N˜m,pq0,kl+∇Δ_εm,pqkl,

_λIFkΔ_φm,pqk−_λIFlΔ_φm,pql=∇Δ_ρpqkl+∇Δ_Tpqkl+_λIF0∇Δ_N˜IF,pq0,kl+∇Δ_εIF,pqkl,

where∇Δ_N˜m,pq0,kland∇Δ_N˜IF,pq0,klare:

∇Δ_N˜m,pq0,kl=2848+^κk2848∇Δ_Nm,pq0,kl−^κl−^κk2848Δ_Nm,pq0,l,

∇Δ_N˜IF,pq0,kl=2848+^κk2848∇Δ_NIF,pq0,kl−^κl−^κk2848Δ_NIF,pq0,l.

The IF combination of∇Δ_NIF,pq0,klandΔ_NIF,pq0,lhas the same formula with (4). As^κl−^κk2848is fairly small (from −0.0046 to +0.0046),Δ_Nm,pq0,landΔ_NIF,pq0,lcould be estimated in advance and do not have to be fixed. In this case, afterΔ_Nm,pq0,landΔ_NIF,pq0,lare applicable, the GLONASS DD ambiguities of (10) and (11) can be estimated and fixed like GPS/Galileo/BDS satellites.

2.3. 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.

2.4. 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:

{_φMW=_φ1−_φ2−_f1−_f2f1+_f2(_P1+_P2)_φGF=_φ1−_λ2λ1φ2.

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. 2.5. 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:

{∇Δ_φMW=∇Δ_N2−∇Δ_N1∇Δ_φGF=1g∇Δ_N2−∇Δ_N1,

whereg=_f2f1. Thus, the DD dual-frequency ambiguity can be calculated by:

{∇Δ_N^1=1g−1∇Δ_φMW+gg−1∇Δ_φGF∇Δ_N^2=gg−1∇Δ_φMW+gg−1∇Δ_φGF.

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:

{∇Δ_N˜1=round(∇Δ_N^1)∇Δ_N˜2=round(∇Δ_N^2),

where∇Δ_N˜1and∇Δ_N˜2are the initial integer ambiguities. Then, the residuals of MW and GF combinations,_rMWand_rGF, after moving the initial integer ambiguities from (13) will be calculated as:

{_rMW=∇Δ_φMW−(∇Δ_N˜2−∇Δ_N˜1)_rGF=∇Δ_φGF−(1g∇Δ_N˜2−∇Δ_N˜1).

Finally, the RMS statistics of the MW and GF combinations of the DD ambiguity arc,_σMWand_σ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). 2.6. 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:

v=A_xa+B_xb+C_xu−l,P,

where the state vector_xacontains coordinate parameters of rover stations;_xbis the vector for ambiguity parameters and vector_xuincludes the other parameters.vdenotes the noise vector.lis the OMC (observed minus computed observations).Pis the weight matrix of DD observations.A,B, andCdenote the corresponding coefficient matrices of the parameters, respectively. Thus, the estimated states and their coefficient matrices can be written into vectorsxandHas:

x=^[_xaxbxu]TH=[ABC].

The parameter states and their covariance matrix will be:

{_xe+1−=_{Φe+1,e xe+Pe+1−}=_{Φe+1,e Pe+ Φe+1,eT}+_Qe,

where_xe+is the a posteriori state vector at epoch e,_xe+1−is the a priori state vector at epoch e+1,_Φe+1,eis the state transition matrix from epoch e to e+1,_Pe+is the a posteriori covariance matrix of_xe+,_Pe+1−is the a priori covariance matrix of_xe+1−, and_Qeis the process noise covariance matrix. The measurement update step is given as:

{_Ke+1=_{Pe+1− He+1T} ^{(_{He+1 Pe+1− He+1T}+_Re+1)−1}_xe+1+=_xe+1−+_Ke+1(_le+1−_{He+1 xe+1−})_Pe+1+=(_Ie+1−_{Ke+1 He+1})_Pe+1−,

where_Ke+1is the Kalman gain matrix,_xe+1+is the a posteriori parameter state vector at epoch e+1,_le+1is the measurement residual vector,_He+1is the measurement design matrix,_Re+1is the measurement noise covariance matrix, and_Ie+1is 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,L1and∇Δ_N^f,L2with the standard deviations_σL1and_σL2in 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 smallest ones. If the float ambiguities and their standard deviations satisfy (21), we can round the float ambiguities as to the fixed ones:

{ratio=^{(1−|dL1|dL1)2}>relrankratio=^{(1−|dL2|dL2)2}>relrank_σL1<_σtσL2<_σt,

where:

{dL1=∇Δ_N^f,L1−round(∇Δ_N^f,L1)dL2=∇Δ_N^f,L2−round(∇Δ_N^f,L2),

whererelrankdefaultly with 25 and_σtwith 0.25 cycles are the thresholds.

In long baselines, the estimated float IF ambiguities and standard deviations are∇Δ_N^f,IFand_σ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_χminto_χ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):

{RL1∈[_χmin,_χmax]RL2∈[RL1−η,RL1+η],

{η=round(|_σMW|+|_γMW|)+1_χmax=round(∇Δ_N^f,LC−βη)/(α+β)+1_χmin=round(∇Δ_N^f,LC+βη)/(α+β)−1.

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 coefficientsaandbin (27) are used to control the effects of MW and GF combinations in the ambiguity resolution. We can adjustaandbto 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:

{dIF=∇Δ_N^f,IF−(αRL1+βRL2)dMW=_rMW−(RL1−RL2)dGF=_rGF−(RL1−_f1f2RL2),

{fi_tIF=^(dIF/_σIF)2fi_tMW=^(dMW/_σMW)2fi_tGF=^(dGF/_σGF)2,

fit=fi_tIF+a·fi_tMW+b·fi_tGF,

ratio=^{(next_fitbest_fit)2}.

If the candidates satisfy all the following conditions, the best candidates will be the fixed ambiguities we searched for. These conditions are:_σIFand_σMWlower 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ξ:

ξ=20/(1−^e−iter5),

whereiteris 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⁹ 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. 3. 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. 3.1. 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.

3.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.

3.2.1. Ambiguity Fixing Performance

Figure 2 and Figure 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.

3.2.2. 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.

3.3. 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.

3.3.1. 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. For the baseline length of about 20 km, the results from the L1+L2 observations and IF combinations have similar RMS statistics, with 5–10 mm in the N and E directions and 15–20 mm in the U direction. For the baselines of NNOR-PERT and GAMG-DAEJ, the precisions of the two baselines can achieve to 6–10 mm in the N and E directions, and 30–40 mm in the U direction. The large residual errors in the U directions for the long baselines are probably because of the residual tropospheric delays.

3.3.2. Residual Analysis

Figure 5 and Figure 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%.

3.4. 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.

3.4.1. 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.

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].

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.

3.4.2. 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.

4. 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.

Figure 1. Data processing flowchart.

Figure 2. Ambiguity resolution of each arc with L1+L2 observations. The experimental data is from PERT-CUT0 on DOY 84 2019. The top panel shows the standard deviations of the float ambiguity estimates of L1 and L2. The middle panel indicates the fractional part of the float ambiguity estimates, and the bottom one presents the lower ratio value of L1 and L2 estimates in (24). The red dash line indicates the threshold. The diamond symbol denotes the length of the arc referring to the right vertical axis. The gray background means the unfixed ambiguity arcs.

Figure 3. Ambiguity resolution of each arc with IF observations. The experimental data is from NNOR-PERT on DOY 84 2019. The top panel shows the standard deviations of the IF estimates referring to the left axis and the MW combinations to the right. The middle panel denotes the fit values computed with (26) and (27). Notes the axis of GF fit values is shown on the right. In the bottom, the ratio values calculated with (28) are shown. The red dash line indicates the lowest threshold of 20. The diamond symbols denote the length of the arc referring to the right vertical axis. The gray background means the unfixed ambiguity arcs.

Figure 4. Residual time series of baseline resolutions after removing averages. Notes the scale ranges of the panels are different.

Figure 5. The posteriori residual RMS and the observation number statistics of GNSS satellites for PERT-CUT0 after solving with L1+L2 observations on DOY 88.

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. Residual distribution histograms of PERT-CUT0 after solving with the L1+L2 observations and IF combinations on DOY 88. The top and middle panels denote the residuals of L1 and L2. The bottom panel shows the IF residuals. The black curve denotes a fitted normal distribution, and⊂means 'within'.

Figure 8. The landscape of the FRB, antenna set-ups with baseline lengths and the observation environment of the SHM2 and SHM4 stations.

Figure 9. Displacement time series of SHM2 and SHM4.

Figure 10. The spectrum analysis of vibrations of SHM2 and SHM4. The purple vertical dash lines indicate the locations of peak frequencies.

Figure 11. Distribution histogram of residuals for SHM2-SHM1.

Figure 12. Trajectory of the UAV after data processing.

Figure 13. Distribution histograms of the residuals after data processing for the UAV flying trial.

Baselines	Locations of Stations	Length	Statement	Data Length, Sampling Rate
CUTA-CUT0	Curtin, Australia	8.42 m	Static	7 d, 30 s
KIRU-KIR0	Kiruna, Sweden	4.47 km	Static	7 d, 30 s
PERT-CUT0	Perth, Australia	22.41 km	Static	7 d, 30 s
NNOR-PERT	New Norcia, Australia	88.56 km	Static	7 d, 30 s
GAMG-DAEJ	Daejeon, Republic of Korea	102.39 km	Static	7 d, 30 s
SHM2-SHM1	Edinburgh, UK	1528.09 m	Kinematic	1 h, 10 Hz
SHM4-SHM1	Edinburgh, UK	1032.67 m	Kinematic	1 h, 10 Hz
ROVE-RS01	Wuhan, China	-	Kinematic	16 m, 5 Hz

Models or Parameters	Strategies
Observations	Short baselines: L1+L2; Long baselines: IF combinations. The dual-frequency observations for different systems, GPS: L1/L2, BDS: B1/B2, GLONASS: G1/G2, Galileo: E1/E5a
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
Cutoff elevation	10°
A priori tropospheric zenith dry/wet delays	Provided by Saastamoinen model [25]. The meteorological parameters are provided by GPT2w model [26]. The mapping function is GMF [27].
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 m2/s.
PCO and PCV	GNSS satellites and receiver antennas follow igs14_2029.atx
Weighting scheme	Elevation dependent model withσ2=a2+b2/sin2θ, whereθpresents the elevation of satellite.
Ephemeris	BRDM combined broadcast ephemeris.

Baselines	Strategy	Sessions (DOY)
		83	84	85	86	87	88	89
CUTA-CUT0	L1+L2	100.0	100.0	100.0	100.0	100.0	100.0	100.0
KIRU-KIR0	L1+L2	100.0	100.0	100.0	100.0	100.0	100.0	100.0
PERT-CUT0	L1+L2	93.24	95.83	96.32	96.96	97.88	95.95	96.18
PERT-CUT0	IF	90.00	93.38	94.53	95.76	96.53	93.65	94.40
NNOR-PERT	IF	91.93	93.33	92.52	91.40	92.34	92.89	90.96
GAMG-DAEJ	IF	92.16	94.48	90.41	92.54	91.43	94.33	91.58
SHM2-SHM1	L1+L2	100.00 (114)
SHM4-SHM1	L1+L2	96.20 (114)
ROV1-RS01	L1+L2	71.76 (88)

Baselines and Directions		Session (DOY)
		83	84	85	86	87	88	89
CUTA-CUT0 (L1+L2, 8.42 m)	N	1.2	1.1	1.2	1.4	1.3	1.3	1.3
	E	1.3	1.2	1.2	1.2	1.3	1.4	1.3
	U	3.0	2.8	2.9	3.2	3.0	3.3	3.0
KIRU-KIR0 (L1+L2, 4.47 km)	N	3.4	3.5	2.6	3.9	3.7	3.9	2.9
	E	3.0	2.6	2.2	3.4	3.5	3.4	3.0
	U	7.0	9.3	7.4	8.4	9.2	9.4	8.6
PERT-CUT0 (L1+L2, 22.41 km)	N	5.5	7.7	8.3	7.8	10.3	9.0	6.8
	E	5.4	5.2	5.1	6.4	5.0	5.1	5.2
	U	15.4	16.1	20.4	19.7	20.8	12.7	16.5
PERT-CUT0 (IF, 22.41 km)	N	5.4	5.2	6.7	6.0	6.9	5.0	7.8
	E	5.9	5.2	5.6	5.7	4.7	5.0	5.6
	U	17.5	14.6	17.6	21.3	22.0	15.3	22.4
NNOR-PERT (IF, 88.56 km)	N	6.6	7.1	8.2	7.0	7.6	8.7	9.7
	E	8.0	8.5	7.9	11.3	7.5	7.2	8.8
	U	34.8	32.2	43.1	39.2	31.8	41.3	43.9
GAMG-DAEJ (IF, 102.39 km)	N	8.3	6.4	6.9	8.7	6.8	9.3	7.1
	E	8.0	6.5	7.9	9.9	8.0	6.9	8.8
	U	28.2	25.5	28.3	32.9	32.3	39.0	33.5

Baselines	Observation	GPS	BDS	GLONASS	Galileo
CUTA-CUT0	L1	4.1	3.9	5.8	4.1
CUTA-CUT0	L2	4.6	4.0	5.2	4.1
KIRU-KIR0	L1	7.5	7.5	8.6	7.9
KIRU-KIR0	L2	10.0	9.0	11.0	9.8
PERT-CUT0	L1	9.2	7.0	12.4	10.3
PERT-CUT0	L2	12.4	9.5	14.2	13.2
PERT-CUT0	IF	12.8	10.1	14.8	13.9
NNOR-PERT	IF	14.9	13.0	17.3	17.0
GAMG-DAEJ	IF	12.6	13.3	14.7	13.4

Stations	N	E	U
SHM2	0.15	0.075 0.268 0.342	0.104 0.205 0.268
SHM4	0.104 0.15 0.18	0.18	0.18

Author Contributions

R.X. proposed the idea and developed the program, R.X. and Q.C. processed the data, wrote this paper and analyzed the results of the experiments; X.M. led the GeoSHM system development; X.M. and W.J. revised the manuscript; X.A. and Q.H. carried out the experiment and programming. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported in part by the National key R&D Program of China under Grant 2018YFC1503601, and the CRSRI Open Research Program (Program SN: CKWV2019749/KY). This research was also supported by the National Natural Science Foundation of China (Grant Nos. 41525014), the National Science Innovation Group Foundation of China (Grant No. 41721003) and the Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University (Grant No. 18-02-05). The European Space Agency also funded the GeoSHM project. Grant numbers 4000108996/13/N4/US and 4000116646/16/NL/US.

Acknowledgments

Thank goes to The Office of China Postdoctoral Council for providing the first author a scholarship which allows him to visit the University of Nottingham for two years to research in the UK from October 2019 (No: 20190050). ESA is acknowledged for sponsoring GeoSHM feasibility a Demo project. Ubipos UK Ltd. is the leading company for GeoSHM Demo Project.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Benassi, F.; Dall'Asta, E.; Diotri, F.; Forlani, G.; di Cella, U.M.; Roncella, R.; Santise, M. Testing accuracy and repeatability of UAV blocks oriented with GNSS-supported aerial triangulation. Remote Sens. 2017, 9, 172.

2. Forlani, G.; Dall'Asta, E.; Diotri, F.; di Cella, U.M.; Roncella, R.; Santise, M. Quality assessment of DSMs produced from UAV flights georeferenced with on-board RTK positioning. Remote Sens. 2018, 10, 311.

3. Li, W.; Li, W.; Cui, X.; Zhao, S.; Lu, M. A tightly coupled RTK/INS algorithm with ambiguity resolution in the position domain for ground vehicles in harsh urban environments. Sensors 2018, 18, 2160.

4. Meng, X.; Roberts, G.W.; Dodson, A.H.; Cosser, E.; Barnes, J.; Rizos, C. Impact of GPS satellite and pseudolite geometry on structural deformation monitoring: Analytical and empirical studies. J. Geod. 2004, 77, 809-822.

5. Meng, X.; Dodson, A.; Roberts, G.W. Detecting bridge dynamics with GPS and triaxial accelerometers. Eng. Struct. 2007, 29, 3178-3184.

6. Meng, X.; Nguyen, D.; Xie, Y.; Owen, J.; Psimoulis, P.; Ince, S.; Chen, Q.; Ye, J.; Bhatia, P. Design and implementation of a new system for large bridge monitoring-GeoSHM. Sensors 2018, 18, 775.

7. Meng, X.; Xi, R.; Xie, Y. Dynamic characteristic of the forth road bridge estimated with GeoSHM. J. Glob. Position. Syst. 2018, 16, 4.

8. Moschas, F.; Stiros, S. Dynamic multipath in structural bridge monitoring: An experimental approach. GPS Solut. 2014, 18, 209-218.

9. Imparato, D.; El-Mowafy, A.; Rizos, C.; Wang, J. Vulnerabilities in SBAS and RTK positioning in intelligent transport systems: An overview. In Proceedings of the International Global Navigation Satellite System Association IGNSS Symposium, Sydney, Australia, 7-9 February 2018.

10. Fan, P.; Li, W.; Cui, X.; Lu, M. Precise and robust RTK-GNSS positioning in urban environments with Dual-antenna configuration. Sensors 2019, 19, 3586.

11. El-Mowafy, A. Analysis of web-based GNSS post-processing services for static and kinematic positioning using short data spans. Surv. Rev. 2011, 43, 535-549.

12. Chen, G. GPS Kinematic Positioning for the Airborne Laser Altimetry at Long Valley, California. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, USA, 1998.

13. Lu, B.; Jin, J.; Duan, W.; Chen, L.; Guan, H. Research of GPS signal multipath effects based on GAMIT TRACK. Adv. Mater. Res. 2012, 588, 912-919.

14. King, M.; Williams, S. Apparent stability of GPS monumentation from short-baseline time series. J. Geophys. Res. 2009, 114.

15. Moschas, F.; Stiros, S. PLL bandwidth and noise in 100 Hz GPS measurements. GPS Solut. 2015, 19, 173-185.

16. Tiryakioglu, I.; Yigit, C.O.; Yavasoglu, H.; Saka, M.H.; Alkan, R.M. The determination of interseismic, coseismic and postseismic deformations caused by the Gökçeada-Samothraki earthquake (2014, Mw: 6.9) based on GNSS data. J. Afr. Earth Sci. 2017, 133, 86-94.

17. Zrinjski, M.; Barković, Đ.; Matika, K. Razvoj i modernizacija GNSS-a. Geodetski list. 2019, 73, 45-65.

18. Teunissen, P.J.G.; Montenbruck, O. Springer Handbook of Global Navigation Satellite Systems; Springer: Berlin, Germany, 2017.

19. Teunissen, P.J.G. A new GLONASS FDMA model. GPS Solut. 2019, 23, 100.

20. Leick, A. GLONASS satellite surveying. J. Surv. Eng. 1998, 124, 91-99.

21. Sleewaegen, J.; Simsky, A.; Wilde, W.; Boon, F.; Willems, T. Demystifying GLONASS inter-frequency carrier phase biases. Inside GNSS 2012, 7, 57-61.

22. Geng, J.; Zhao, Q.; Shi, C.; Liu, J. A review on the inter-frequency biases of GLONASS carrier-phase data. J. Geod. 2017, 91, 329-340.

23. Xi, R.; Meng, X.; Jiang, W.; An, X.; Chen, Q. GPS/GLONASS carrier phase elevation-dependent stochastic modelling estimation and its application in bridge monitoring. Adv. Space Res. 2018, 62, 2566-2585.

24. Gurtner, W.; Estey, L. RINEX: The Receiver Independent Exchange Format, Version 3.02; Technical Report; IGS Central Bureau: Pasadena, CA, USA, 2013.

25. Saastamoinen, J. Atmospheric correction for the troposphere and stratosphere in radio ranging satellites. Use Artif. Satell. Geod. 1972, 15, 247-251.

26. Böhm, J.; Möller, G.; Schindelegger, M.; Pain, G.; Weber, R. Development of an improved empirical model for slant delays in the troposphere (GPT2w). GPS Solut. 2015, 19, 433-441.

27. Böhm, J.; Niell, A.; Tregoning, P.; Schuh, H. Global Mapping Function (GMF): A new empirical mapping function based on numerical weather model data. Geophy. Res. Lett. 2006, 33, L07304.

28. King, M.; Watson, C.; Penna, N.; Clarke, P. Subdaily signals in GPS observations and their effect at semiannual and annual periods. Geophys. Res. Lett. 2008, 35, L03302.

29. King, M.; Watson, C. Long GPS coordinate time series: Multipath and geometry effects. J. Geophys. Res. Sol. Earth 2010, 115, B4.

30. Xi, R.; Zhou, X.; Jiang, W.; Chen, Q. Simultaneous estimation of dam displacements and reservoir level variation from GPS measurements. Measurement 2018, 122, 247-256.

31. Meng, X.; Nguyen, D.; Owen, J.; Xie, Y.; Psimoulis, P.; Ye, G. Application of GeoSHM System in Monitoring Extreme Wind Events at the Forth Road Bridge. Remote Sens. 2019, 11, 2799.