Integrity Monitoring for Horizontal RTK Positioning: New Weighting Model and Overbounding CDF in Open-Sky and Suburban Scenarios

: Integrity monitoring is an essential task for ensuring the safety of positioning services. Under a selected probability of hazardous misleading information, the protection levels (PLs) are computed according to a considered threat model to bound the positioning errors. A warning message is sent to users when the PL exceeds a pre-set alert limit (AL). In the short-baseline real-time relative kinematic positioning, the spatially correlated errors, such as the the orbital errors and the atmospheric delays are signiﬁcantly reduced. However, the remaining atmospheric residuals and the multipath that are not considered in the observation model could directly bias the positioning results. In this contribution, these biases are analysed with the focus put on the multipath effects in different measurement environments. A new observation weighting model considering both the elevation angle and the signal-to-noise ratios is developed and their impacts on the positional results are investigated. The coefﬁcients of the proposed weighting model are determined for the open-sky and the suburban scenarios with the positional beneﬁts maximised. Next, the overbounding excess-mass cumulative distribution functions (EMCs) are searched on the between-receiver level for the weighted phase and code observations in these two scenarios. Based on the mean and standard deviations of these EMCs, horizontal protection levels (HPLs) are computed for the ambiguity-ﬁxed solutions of real experiments. The HPLs are compared with the horizontal positioning errors (HPEs) and the horizontal ALs (HALs). Using the sequential exclusion algorithm developed for the ambiguity resolution in this contribution, the full ambiguity resolution can be achieved in around 100% and 95% of the time for the open-sky and the suburban scenarios, respectively. The corresponding HPLs of the ambiguity-ﬁxed solutions are at the sub-dm to dm-level for both scenarios, and all the valid ambiguity-ﬁxed HPLs are below a HAL of 0.5 m. For the suburban scenario with more complicated multipath environments, the HPLs increase by considering extra biases to account for multipath under a certain elevation threshold. In complicated multipath environments, when this elevation threshold is set to 30 degrees, the availability of the ambiguity-ﬁxed solutions could decrease to below 50% for applications requiring HAL as low as 0.1 m.


Introduction
The real-time kinematic (RTK) positioning technique is a popular Global Navigation Satellite System (GNSS) based positioning technique.It is substantially different from the precise point positioning (PPP) technique that typically requires a long convergence time of tens of minutes [1][2][3], because the RTK positioning technique removes or mitigates observation errors such as clocks, hardware biases, and a large part of the spatially correlated errors such as the atmospheric delays.This enables ambiguity resolution within a short time.For the RTK of short baselines, e.g., less than 10 km, the atmospheric delays can be considered significantly reduced under quiet atmosphere conditions.The remaining biases mainly consist of the multipath effects.They could vary strongly with the measurement environments, for instance in urban areas, and thus may hamper successful ambiguity resolution and bias the positioning results.
To de-weight the observations with large noise, multipath and atmospheric residuals, elevation-dependent weighting functions are typically used in the RTK positioning [4].The elevation angle, however, is not the only factor that could be used to build a weighting model.As shown in [5,6], the weighting model exploiting the carrier-to-noise density ratio (C/N 0 ) can also effectively de-weight low-quality observations and improve the positioning results.Kuusniemi et al. [7] also utilised C/N 0 to weight the code measurements in signal-degraded environments.In environments with complicated multipath effects and non-line-of-sight (NLOS) errors, e.g., in urban canyons, considering both the elevation angles and the C/N 0 in the weighting model was found to be useful for code-based positioning [8].Medina et al. [9] also used mixed elevation-based and C/N 0 -based model to weight the code noise and multipath in signal-degraded scenarios.Further studies were also performed with respect to the stochastic modelling of baseline data [10][11][12].
In RTK positioning, the multipath effects remaining in the double-differenced (DD) observations strongly vary with the measurement environments.This variation also differs for different observation types.For a weighting model considering both the elevation angles and the C/N 0 , it is thus important to adjust the impacts of these two factors under different measurement scenarios and for different observation types.In this contribution, the first issue to be addressed is to develop a weighting model allowing different and adjustable impacts of the elevation angles and the C/N 0 .The weighting coefficients of these two factors are searched for the open-sky and the suburban scenarios, under the ambiguity-float and -fixed modes.For the ambiguity-fixed case, a sequential exclusion algorithm is developed for ambiguity resolution, so that possibly many ambiguities at different elevation angles and C/N 0 can be resolved at each epoch, and the remaining residuals can be used for the determination of the weighting coefficients.The weighting models achieving the largest positional benefits are then used for further analysis.
In road transport applications, integrity monitoring (IM) is an essential task to ensure that the positioning errors are less than the alert limit (AL) with a pre-defined probability of hazardous misleading information (PHMI).With the varying satellite geometry over time, the protection levels (PLs) are computed to bound possible positioning errors under the PHMI, and are compared with the AL.In case that the AL is exceeded, a warning message is sent to users within the time-to-alert.The IM was investigated for aviation over several decades [13][14][15][16][17].In one IM approach, the solution separation test (SST) for the fault detection and exclusion (FDE) procedure is first applied.The protection levels are produced next within the advanced receiver autonomous integrity monitoring (ARAIM) algorithm for both the horizontal and vertical guidance [13,14,18].The ARAIM algorithm in aviation has been adapted in recent years for ground-based applications such as road transport.The IM was performed, e.g., for PPP [19] and RTK [16,20,21].
The ARAIM algorithm used in aviation, however, is applied for single-receiver positioning.For the ground-based RTK positioning, more complicated multipath environments may need to be considered in the threat model, while the accuracy of the satellite clocks and orbits, which are important for the aviation ARAIM, do not play an important role for short-baseline RTK processing.In addition, the atmospheric residuals need to be considered at the between-receiver level.In this contribution, the analysis is performed for short baselines in atmosphere quiet days, so that the focus is on the multipath effects under different measurement scenarios.This leads to the second issue to be addressed in this contribution, i.e., determination of the overbounding mean and standard deviations for the between-receiver phase and code noise/biases under different measurement scenarios, which will be used for computation of the representative PLs during the IM.It is noted that during geophysical events like solar flare [22] and magnetic storms [23,24], abnormal ionosphere behaviours may lead to enlarged differential ionospheric errors in the RTK processing.Note that this study is focused on short-baseline RTK positioning in quiet atmosphere conditions, and does not attempt to account for the RTK integrity monitoring under extreme atmospheric conditions.
To solve this problem, the excess-mass cumulative distribution functions (EMCs) are searched for different observation types and measurement scenarios.In this contribution, the authors do not address the very challenging urban areas, which will be considered in future research, but put the focus on open-sky and suburban scenarios.In the presented IM approach, after applying the FDE procedure with the SST, utilising these overbounding mean and standard deviations, the horizontal protection levels (HPLs) are computed for the ambiguity-fixed solutions.Using the sequential exclusion method for ambiguity resolution as mentioned before, the proportion of the time epochs achieving full ambiguity resolution (FAR) is discussed for the open-sky and the suburban scenarios.The valid ambiguity-fixed HPLs are checked to bound the horizontal positioning errors (HPEs) and to be less than the horizontal AL (HAL).In addition, for the suburban scenario having more complicated multipath environments, investigations are also performed considering extra multipath biases for low-elevation rover phase observations.The resultant HPLs are compared with a set of different HALs that might be needed for different road transport applications.
The paper starts with the signal analysis of different observation types under open-sky and suburban scenarios.The weighting model considering both the elevation angles and the C/N 0 is developed, and the coefficients leading to the largest positional benefits are given and used in the positioning and IM.The overbounding EMCs are then searched for the between-receiver mode for different observation types and various scenarios.Subsequently, the threat model, the FDE method, and the procedure to compute the HPLs are explained in detail.Test results based on real data are then analysed, where the positioning solutions, the HPLs and the HALs are discussed.The conclusions are then drawn.

Signal Analysis
For the GPS dual-frequency scenario, the expected observed-minus-computed (O-C) vectors of the DD short-baseline phase (∆Φ DD ) and code (∆P DD ) observations can be formulated as [25]: where to the rover receiver, and the differencing operator D T m = [−e m−1 , I m−1 ] forms the between-satellite difference with e m−1 and I m−1 denoting the vector of ones with the length m − 1 and the identity matrix with size (m − 1) × (m − 1), respectively.m denotes the number of satellites with 1 representing here the reference satellite.λ j stands for the wavelength on frequency j, and diag(•) denotes the diagonal matrix with the diagonal elements contained in (•).E(•) is the expectation operator.The rover coordinate increment vector starting from an approximate position and the DD ambiguity vector are represented by ∆r and a, respectively.The operator ⊗ denotes the Kronecker product.Apart from the a priori DD receiver-satellite geometric distance, the corrections of the antenna phase centre offset (PCO) and phase centre variation (PCV), and the a priori DD tropospheric delays computed based on the Saastamoinen model [26] and the Ifadis mapping function [27] are also corrected for the O-C terms.
Assuming that the code and phase observations are uncorrelated, the variance-covariance matrix of the O-C terms, Q y can be formulated as: where the variance-covariance matrices of the between-receiver phase (Q φ ) and code (Q p ) observations are diagonal matrices, given as: for which σ φ,j and σ p,j represent the between-receiver phase and code standard deviations on frequency j in the zenith direction (and for a between-receiver C/N 0 of ∞), respectively.blkdiag(•) forms the block diagonal matrix.The phase and code weighting matrices W φ,j and W p,j are constructed as follows: for which the phase and code weighting functions w s φ,j and w s p,j for satellite s on frequency j will be discussed later.
For the purpose of signal analysis, and with the ground truth of the baseline coordinates known, the DD phase (δφ 1s DDj ) and code (δp 1s DDj ) residuals can then be expressed as: where ∆φ 1s DDj and ∆p 1s DDj represent the DD O-C terms of the phase and code observations on frequency j, respectively, and ∆ρ 1s DD denotes the DD geometry increment based on the ground truth and a priori coordinates of the rover.Using the strong baseline-known model [28], the DD ambiguities N 1s j are resolved with the search-and-shrink least-squares ambiguity decorrelation adjustment (LAMBDA) method using the integer least-squares (ILS) estimator [29][30][31].The ambiguities are only considered as successfully resolved when the ratio µ = R 2 /R 1 is higher than 3 [32,33].The R 1 and R 2 represent here the squared norm of the ambiguity residuals measured in the metric of the variance-covariance matrix of the float ambiguities for the best and the second best integer ambiguity candidates, respectively.The resolved ambiguities are assumed fixed until cycle slips or large data gaps are detected, and the unresolved ambiguities are treated as pseudo-observations and constrained as constant in time.To resolve possible large ambiguity sets under different elevation angles and C/N 0 , so that possibly many phase residuals (Equation ( 7)) can be utilised in the data analysis, a sequential exclusion algorithm is developed and applied to exclude suspected observations that may hamper the ambiguity resolution.The detailed procedure is described in Appendix A. With the baseline coordinates and the ambiguities corrected in Equations ( 7) and ( 8), the DD residuals should contain the noise, the multipath, and the DD atmospheric residuals.
The data analysis is performed using GPS dual-frequency data with 1 Hz sampling rate collected on L1C and L2W channels.To simulate different measurement scenarios, the continuously operating reference station (CORS) UWA0 located at the University of Western Australia is used as the reference station, and the station CUT0 located on the roof of a building in the Curtin University (left panel of Figure 1), which is about 8 km away from UWA0, is used as the rover station to simulate the open-sky scenario.With similar baseline length, another rover receiver mounted on the roof of a lower building, which is surrounded by trees and metal structures (middle and right panels of Figure 1), is used to simulate the suburban scenario.The station is denoted as SUB0.The receiver and antenna types of both baselines are given in Table 1.The data collected from 6:34:45 to 23:59:59 in GPS time (GPST) on 15 October 2019 were analysed for both baselines.The data were collected under quiet ionosphere conditions.Based on the Australian total electron content (TEC) maps provided by the Australian Bureau of Meteorology/Space Weather Branch (IPS) [34], during the test time periods, the test region had TEC values lower than 20 TEC unit (TECU).For the troposphere, the daily rainfall amount is 0 mm [35].In this study, the lower bound of the C/N 0 in the contribution is set to 12 dBHz, which is recognised as the minimal possible signal strength [36].The satellite having the maximum mean elevation angles for the baseline is used as the reference satellite.It remains the reference satellite until its mean elevation angle is reduced to below 30 degrees.The elevation mask is set to 10 degrees in this study.It is noted that the reference station of RTK processing is assumed to be located under good measurement environment and broadcasts its data in a continuous manner.

Elevation, C/N 0 and Biases
To show the changes of the DD residuals with elevation angles, the elevation angle θ 1s is computed based on the elevation-dependent weighting function w(θ s ) = 1/ 1 + 10 × exp(−θ s /10 • ) 2 (see also Appendix A), so that: where σ j,0 represents the between-receiver zenith-referenced standard deviation for code or phase observations on frequency j, and σ 1s j represents that on the DD level, mapped to the slant direction.Based on Equation ( 9): so that θ 1s can be obtained with: where log(•) is the natural logarithm of the element in (•).
The DD phase (δφ 1s DDj ) and code (δp 1s DDj ) residuals are grouped for the elevation intervals of θ 1s for every one degree, and the root mean square (RMS) of the residuals are plotted in the top panels of Figure 2 for the open-sky and suburban scenarios.The code RMS is divided by a factor of 100 for the purpose of illustration.Together with the RMS, the total number of observations used for computing the RMS in each elevation interval is also shown in the bottom panels of Figure 2. In general, the RMS of all observation types decrease with increasing elevation angles.Compared to the open-sky scenario, larger increase in the phase RMS can be observed for the suburban scenario when decreasing the low elevation angles (the red and blue lines in the top panel of Figure 2b).This should be related to the larger phase multipath effects in the suburban scenario.For the code RMS, compared with their smooth behaviours in the open-sky scenario, greater randomness in RMS is visible under the suburban scenario due to the fact that the residuals become less sensitive to the elevation angles in the more complicated environment (the green and magenta lines in the top panel of Figure 2b).In addition to the relationship between the DD residuals and the elevation angles, the changes of the residuals with the C/N 0 were also analysed.Similar to Equation ( 9), the C/N 0 S 1s j are computed based on the C/N 0 -dependent weighting function w(S s j ) = 1/(1 + b j × 10 −0.1×S s j ) (see also [6]), such that: where the coefficient b j varies with the frequency, and S s j denotes the between-receiver C/N 0 for satellite s on frequency j.Based on Equation ( 12): so that S 1s j can be obtained with: where the between-receiver C/N 0 S s j can be obtained, similar to Equation ( 14), with: where S s j,r and S s j,v represent the C/N 0 of the reference and the rover receivers for satellite s on frequency j, respectively.
The RMS of the DD residuals is computed within small C/N 0 intervals (S 1s j ) of 0.3 dBHz starting from the minimum S 1s j above the lower bound of 12 dBHz, and is shown in the top panels of Figure 3.In general, a higher correlation between the RMS and the C/N 0 can be observed in the suburban scenario.It can also be observed that the signals on L2W (the red and magenta lines in the top panels) generally experience lower C/N 0 than those on L1C (the blue and green lines in the top panels).This agrees with the results shown, e.g., in [37] that the C/N 0 for the L2 P(Y) signals are lower than those for the L1 C/A and L5 signals, which is related to the lower minimum received power of L2.
Based on the analysis above, the weighting functions w s φ,j and w s p,j are in this study considered as a function of both the elevation angles θ s and the C/N 0 for phase S s j in the form: where the parameters a φ and a p are used to adjust the elevation impact on the weights, and b φ,j and b p,j are used to adjust the impact of the C/N 0 on the weighting functions.In this contribution, to evaluate the weighting model as will be discussed in Section 2.2, a φ and a p are empirically searched from 0 to 16 with a step of 2, and b φ,j and b p,j are searched for 0, 100, 500, 1000, 4000, 7000 and 10 4 Hz to find the most appropriate coefficients.This study does not attempt to search for higher coefficients for b φ,j and b p,j , as the weighting functions 1/(1 + b φ,j × 10 −0.1×S s j ) and 1/(1 + b p,j × 10 −0.1×S s j ) are expected to be around one at a S s j of 54 dBHz, which is characterised as the maximum possible signal power strength according to [36].Note that the RMS of the code residuals is divided by a factor of 100.In the bottom panels, the blue and the red lines are almost overwritten by the green and the magenta lines.

The Weighting Function
This section searches for the appropriate coefficients a φ , a p , b φ,j and b p,j in the proposed phase and code weighting functions w s φ,j (Equation ( 16)) and w s p,j (Equation ( 17)).As shown in Figures 2  and 3, the DD residuals are correlated with both the elevation angles and the C/N 0 .In actual GNSS positioning, however, too much de-weighting of the observations at low elevation angles or with low C/N 0 is nearly equivalent to excluding these observations.In such cases, although their large biases will less influence the positional results, the observation model also becomes weaker, especially in the measurement environment with limited satellite visibility.In this section, the authors directly evaluate the influences of the DD residuals on the positioning results.The first case utilizes the dual-frequency phase residuals δφ 1s DDj and approximates the situation in the ambiguity-fixed case, where the code observations are strongly de-weighted and their weighting coefficients thus do not play much of a role in both the open and suburban scenarios with non-extreme code multipath.The second case utilises the dual-frequency code residuals δp 1s DDj and approximates the situation in the ambiguity-float case (at the first epoch of initialisation), where the phase observations are totally reserved for the ambiguities [28] and their weighting coefficients thus do not influence the positioning results.Their influences on the positional results ∆r φ in the first case, and ∆r p in the second case, are calculated from the least-squares adjustment as follows: where the matrix A = e 2 ⊗ (D T m G) (Equation ( 1)), and note that G contains the satellite-to-receiver unit vectors in the north, east and vertical directions.The phase and code residual vectors δφ DD and δp DD consist of the phase and code residuals δφ 1s DDj (Equation ( 7)) and δp 1s DDj (Equation ( 8)) for j = 1, 2. The zenith-referenced undifferenced signal standard deviations are first assumed as 0.003 m and 0.3 m for the phase and code observations, and the between-receiver phase and code standard deviations σ φ,j (Equation (3)) and σ p,j (Equation ( 4)) are obtained by multiplying them by a factor of √ 2 due to the between-receiver differencing.An outlier detection procedure is performed using the χ 2 -test.The test statistics T φ and T p for the phase and code processing are expressed as [13]: The T φ and T p are compared with the inverse cumulative distribution function (CDF) of the central , where f denotes the degrees of freedom, and P FA,χ 2 is the probability of false alert allocated to the χ 2 -test, which is set to 10 −6 in this study.Data with the PDOP higher than 10 is not considered.The PDOP for RTK is generated using the design matrix and the weighting matrix for observations on L1 [28].We note that the χ 2 -test is performed separately for the code-and phase-processing.As the residuals are already corrected by the ambiguities and the ground truth of the rover (Equations ( 7) and ( 8)), they are assumed to contain only the noise and miss-modelled effects like the multipath and the remaining atmospheric residuals.
To evaluate the impacts of the weighting coefficients on the positional results, the weighting coefficients are searched as mentioned at the end of Section 2. The mean of the 3D positioning errors, i.e., the mean of |∆r φ | for phase and |∆r p | for code, is used for the purpose of evaluation.| • | denotes the norm of the vector.The weighting coefficients leading to the smallest mean are given for the phase and code solutions, and for different measurement scenarios in Table 2.The improvement ratio P is computed as: Mean φ (0, 0, 0) , Phase processing Mean p (0, 0, 0) − Mean p ( âp , bp,1 , bp,2 ) Mean p (0, 0, 0) , Code processing (22) where Mean φ (u 1 ,u 2 ,u 3 ) and Mean p (u 1 ,u 2 ,u 3 ) refer to the mean of the phase-and code-based 3D positional errors applying the three weighting coefficients u 1 , u 2 and u 3 , respectively, and âφ , âp , bφ,j and bp,j denote the weighting coefficients leading to the smallest mean value for the corresponding observation types.The case when setting all coefficients to zeros represents the case without using any weighting function.Note that a φ and a p , used as weighting coefficients for the elevation angles, are not frequency-dependent, whereas b φ,j and b p,j , used as weighting coefficients for the C/N 0 , are frequency-dependent.The coefficient combinations delivering invalid solutions, i.e., with PDOP>10 or with too many observations screened in the χ 2 -test, in more than 1% of the tested time epochs are not considered for further evaluations.To make a fair comparison, only the time epochs that produce valid positional results using all valid coefficient combinations are used for computing the statistics in Table 2.For comparison purposes, the statistics are also given for the following cases: • When bφ,j and bp,j (j=1, 2) are set to zero, i.e., using only the elevation-dependent weighting function in the model.

•
When âφ and âp are set to zero, i.e., using only the C/N 0 -dependent weighting function in the model.
Together with the improvement ratio P, the actual improvement in range, denoted as P m , is also given in Table 2 with P m = P × Mean φ (0, 0, 0) for phase processing, and P m = P × Mean p (0, 0, 0) for code processing.
As shown in Table 2, using the proposed weighting model related to both the elevation angle and the C/N 0 , larger positional improvements can be achieved in the suburban scenario than in the open-sky scenario.Higher elevation-related weighting coefficients, i.e., about 10 to 16, are needed for the phase processing than those for the code processing, i.e., about 0 to 6 for the latter.Compared with the elevation-dependent weighting function with b φ,j and b p,j set to zero, the C/N 0 -related weighting functions can generate larger improvements in some cases.For the code processing in the suburban scenario as an example, using solely the C/N 0 -related weighting function has led to an improvement of about 19%, which is much higher than using solely the elevation-dependent weighting function with an improvement of about 8%.16) and ( 17)) leading to the smallest mean value in the positional results and the improvement ratios P of the mean of 3D positioning errors compared to the case that all the weighting coefficients are set to zero (i.e., no weighting model used).P m refers to the actual improvement in the mean of the 3D positional errors.Figure 4 shows the variation in the improvement P m in range with the coefficients a φ for the phase processing and a p for the code processing.The dashed lines represent the case using only the elevation-dependent weighting functions, i.e., setting b φ,j and b p,j to zero, and the solid lines represent the case when setting b φ,j and b p,j to the C/N 0 coefficients leading to the largest P m for the corresponding a φ and a p , respectively.The differences between the solid and the dashed lines for the same scenario, i.e., the same colour, denote thus the further contribution of the C/N 0 to the elevation-dependent weighting models.It can be observed that for both the phase and the code processing, incorporating the C/N 0 in the weighting model could generally lead to larger contributions in P m in the suburban scenario (see the red lines) than in the open-sky scenario.From Figure 4, the highest value of the elevation-related weighting coefficient is observed for the phase processing in the suburban scenario.This is consistent with the information shown in Figure 2 (the blue and red lines in the top panel of Figure 2b).Using the weighting coefficients giving the largest positional improvements (Table 2), the between-receiver phase (σ φ,j ) and code (σ p,j ) standard deviations in the zenith direction and at the S s j of ∞ can be obtained with the least-squares variance component estimation (LS-VCE) procedure [38], and are given in Table 3 for both frequencies and both measurement environments.Table 3. Between-receiver phase (σ φ,j ) and code (σ p,j ) standard deviations at the zenith direction and with the C/N 0 of ∞.

Positioning Method
Before resolving the ambiguities, ambiguity-float RTK positioning is performed using the recursive least-squares algorithm constraining the ambiguities in time [39].Based on the observation model given in Equation ( 1), the time-updated ambiguity vector from the time point t i−1 , denoted as â(t i|i−1 ), and its variance-covariance matrix Q â(t i|i−1 ) can be expressed as: where the matrix F identifies the time-updated ambiguities from the estimated unknown vector ∆ x at t i−1 .Q ∆ x denotes the variance-covariance matrix of ∆ x.No system noise is added in Equation ( 24), as the ambiguities are constrained as constants in time.
The estimates of the unknown vector ∆ x(t i ) can be obtained by adding the estimated ambiguities â(t i|i−1 ) as pseudo-observations in the observation equations at t i and solving using the least-squares adjustment.The updated observation equations and the corresponding variance-covariance matrix are then written as: The ambiguity-float solutions ∆ x at t i can be obtained with the least-squares adjustment as: As mentioned earlier, using the search-and-shrink LAMBDA method with the ILS estimator [29][30][31], the float ambiguities are considered as successfully fixed when the ratio test µ = R 2 /R 1 > 3 (Section 2) passes.The resolved ambiguities are fixed to the obtained integer values ǎ until cycle slips or large data gaps are detected.In cases when the ratio test does not pass, the sequential exclusion algorithm is performed, which is described in Appendix B.
When all ambiguities of the current epoch are fixed, the ambiguity-fixed positional increment ∆ř can be determined as: with (Equation (1)) To demonstrate the effectiveness of the proposed ambiguity resolution method, Figure 5 shows the L1 and L2 residuals for all observed satellites above the elevation mask introducing the resolved ambiguities based on the model presented in this section and the ground truth of the baselines computed from the known positions of the points under both the open-sky and the suburban scenarios.The time period 0:00:00-6:56:01 (GPST) on 17 October 2019 was used for the computations.As shown in Figure 5a, the residuals in the open-sky scenario are generally within ±0.2 cycles, while those in the suburban scenario are larger due to the increased impact of multipath.Benefiting from the sequential exclusion algorithm and the fix-and-hold strategy applied for ambiguity resolution, the time epochs achieving successful FAR has increased from about 80% to 95% in the suburban scenario compared to the case without using the sequential exclusion algorithm.For the time epochs with partial ambiguity resolution (PAR), the average fix-rate, i.e., the averaged number of the fixed ambiguities among all the possible ambiguities, has also increased from about 55% to 91% upon applying the sequential exclusion algorithm.In the open-sky scenario, the FAR was achieved in almost all epochs, i.e., above 99.8%,whether or not the sequential exclusion algorithm is used.Note that the PAR aims to resolve a part of the ambiguities in case that the FAR cannot be achieved.Different from the PAR based on a pre-defined formal ambiguity success rate (ASR) [40,41], the PAR in this study refers to the resolution of a subset of the remaining float ambiguities, which allows the ratio test to be passed.The selection of the subset was determined based on the sequential exclusion algorithm discussed before and given in Appendices A and B for the baseline-known and -unknown cases, respectively.The algorithm is developed to deal with the miss-modelled effects that are difficult to be covered by the model and the formal ASR.

Integrity Monitoring
This section presents the approach for the IM for the ambiguity-fixed RTK positioning under quiet atmosphere conditions.The mean values and standard deviations of the overbounding EMCs are first searched for phase and code residuals on the between-receiver level for both frequencies.Using the determined overbounding EMCs, the threat model protects against the onset of satellite faults, which is assumed to have the probability 10 −5 per satellite per approach [42,43].Assuming that both receivers are tracking dual-frequency phase and code signals from the same m satellites, in addition to the all-in-view solutions, in the DD mode, m − 1 subset solutions can be computed excluding all observations from one satellite in each subset solution.

Overbounding CDF
Due to the existence of multipath and the atmospheric residuals, the empirical distributions of the between-receiver phase and code residuals are neither symmetric nor have a zero-mean.In the IM, as it is difficult to model the exact empirical distribution, an overbounding normal distribution is thus needed, so that the CDF overbounds that of the empirical distribution.In this contribution, the overbounding CDF is searched for the weighted between-receiver phase (δφ v,j ) and code (δp v,j ) residuals, which are formulated as: where the DD weighting functions w 1s φ,j and w 1s p,j are computed applying the variance rule: In such a case, the mean ( mφ,j , mp,j ) and standard deviations ( σφ,j , σp,j ) of the searched overbounding CDF refer to the between-receiver level at the zenith direction and with a C/N 0 of ∞.Note that the residuals screened out during the χ 2 -test (Equations ( 20) and ( 21)) and the time epochs with PDOP higher than 10 are not used for computing the overbounding CDF, as these observations are excluded in the FDE procedure of the RTK processing.
According to [44], the EMC overbounding strategy does not only have no shape restriction on the empirical distribution, but is also evidenced for the conservatism during the transformation from the observation-domain to the position-domain.Compared with the general paired overbounding strategy [45], the EMC allows an excess mass , i.e., a total mass 1 + over 1.This relieves the strict search process at the tails.The empirical distribution, denoted by G a , is then bounded between the left overbounding CDF G L and the right overbounding CDF G R , such that it is always below the G L and above the G R , which are given as: for which the normal distribution with mean m and standard deviation σ is denoted by N ( m, σ).As examples, Figure 6 shows the empirical CDF G a and its overbounding CDFs G L and G R for the weighted code residuals on L1 δp v,1 in the open-sky (a) and the suburban (b) scenarios.To satisfy Equations ( 35) and ( 36), for code residuals, the overbounding mean and standard deviations are searched from 0.05 m to 1 m with a step of 0.01 m, and for phase residuals from 0.001 m to 0.03 m with a step of 0.001 m.The excess mass is set to 0.01 in this contribution and will be considered by computing the protection levels in Section 4.3.The search begins from the smallest overbounding mean value, and for each mean value, the algorithm searches again from the smallest overbounding standard deviation.The search stops when Equations ( 35) and (36) are fulfilled.The overbounding mean and standard deviations for the weighted between-receiver phase and code residuals are given in Table 4 for all observation types, for the open and suburban scenarios.The overbounding mean and standard deviations are in general higher in the suburban scenario than in the open-sky scenario.Table 4. Overbounding mean values ( mφ,j , mp,j ) and standard deviations ( σφ,j , σp,j ) for the weighted between-receiver phase (δφ v,j ) and code (δp v,j ) residuals.The multipath under the open-sky and the suburban scenarios, as well as the atmospheric residuals under quiet conditions are included in the weighted between-receiver residuals (δφ v,j and δp v,j ) and are covered by the overbounding standard deviations and mean values.They will be considered by computation of the protection levels.

FDE Procedure
The FDE procedure used for the RTK is borrowed from the ARAIM algorithm and consists of the SST that is applied in the position domain and a confirmation check (sanity check) with the χ 2 -test applied in the observation domain [13,14,42,46].In this contribution, the positioning integrity monitoring is performed for the ambiguity-fixed solutions having the SST and the χ 2 -test performed before computing the protection levels.The χ 2 -test is also performed for the ambiguity-float solutions to exclude large outliers.
The SST is applied to exclude outliers by comparing the differences between the all-in-view solution and the subset solution excluding suspected satellites with a test threshold.The threshold is calculated based on the formal standard deviation of the corresponding solution difference and the probability of false alert.The ambiguity-fixed all-in-view solution is calculated with ∆ř 0 = S 0 y x as shown in Equation (28), and the corresponding subset solution ∆ř k is computed excluding certain observations (k) using the least-squares adjustment: with where the rows in A x k and the columns and rows in Q y k related to the observations k are excluded during the calculation, and the corresponding column of zeros are then removed in A x k .Note that these columns of zeros for the observations k are added in S k to match its size with that of S 0 .Assuming that the observation errors are normally distributed, and accordingly the positioning errors, the test thresholds of SST in the north (q = 1), east (q = 2) and up (q = 3) directions are expressed as: with where vector o q contains one for the q-th positioning element, and its other elements are set to zero.
N is the number of fault modes, and Q −1 (•) is the inverse right-folded CDF of a standard normal distribution.In contrast to aviation, the focus in ground-based applications is on the horizontal positioning.As a result, in this study, the probability of false alert allocated in the horizontal (P FA,H ) and the vertical (P FA,V ) domains are set to 3 × 10 −6 and 10 −6 , respectively.Note that these values are not yet clearly defined for diverse road transport applications.The SST passes, when for q = 1, 2 and 3.For all the subset solutions that cannot pass the SST, the subset delivering the largest value of |o T q (∆ř 0 − ∆ř k )|/T k,q (for q = 1, 2 or 3) is determined, and the corresponding observations k are excluded.The SST is re-computed until all subsets passing the test is identified.
After the SST, a confirmation test in the observation domain is performed using the χ 2 -test.The test passes, when Recall that the probability of false alert allocated to the χ 2 -test is set to 10 −6 in this study.In case that the χ 2 -test does not pass but the SST passes, the protection levels computed in the next subsection are not considered valid.

Protection Level
Protecting against the onset of satellite faults, the probability of two or more faults is calculated as [14]: where p j represents the probability of the onset of satellite fault for satellite j.N is the number of the fault modes, here the number of the satellites minus 1.Note that as only the GPS observations are used, the constellation fault and fault of the reference satellite are not feasible for a consistency check, and are thus not protected in this study.The P n amounts to, e.g., about 10 −9 to 4.5 × 10 −9 in the case that N equals 5 to 10. Assuming a horizontal probability of hazards misleading information (PHMI H ) of 10 −5 per approach, the P n , which is considered as the probability of non-monitored faults, is at an insignificant level.The probabilities of the monitored faults in the north (q = 1) and east (q = 2) directions, could then be computed as: As in the ARAIM algorithm [13,14,42], the north (q = 1) and east (q = 2) protection levels can be computed (as given in these references): where n 0 and n k are the numbers of the between-receiver observations in the all-in-view and subset cases (k), respectively.bk,q and σk,q , including the all-in-view case, i.e., k = 0, are calculated based on the overbounding standard deviations and mean values as follows (Section 4.1): where diag(•) retrieves the diagonal elements of the matrix contained in (•) as a vector.Recall that the is the excess mass of the EMC (Equations ( 35) and ( 36)), and a total mass 1 + needs to be considered for all the between-receiver observations [44].Qy k is calculated with Equations ( 2)-( 4) excluding the rows and columns for observations k and using σφ,j and σp,j instead of σ φ,j and σ p,j .The vector mDD,k contains the maximum nominal biases for the DD observations, and can be expressed as: mDD,k = mφ,1 , mφ,2 , mp,1 , mp,2 where the items corresponding to observations k are excluded.
Assuming that the PHMI q is equally allocated in each fault mode (including the all-in-view mode), an overbound of the PL q can be computed as the maximum of the PLs computed for each mode [18,19,42,43]: where PL k,q is calculated from PL 0,q = K 0,q × σ0,q + b0,q (52) PL k,q = K k,q × σk,q + bk,q + T k,q , k = 0 (53) for which Finally, the HPL, which bounds the overall horizontal positioning error, is computed as: for which PL 1 and PL 2 are the PLs in the north and east directions, respectively.

Test Results
The baselines UWA0-CUT0 and UWA0-SUB0 (Figure 1) were used for the tests in the open-sky and the simulation of suburban scenarios, respectively, based on 1 Hz GPS data.To test the positioning behaviours under different satellite geometry conditions, two time periods were used for the analysis, i.e., 0:00:00-6:56:01 (GPST) on 17 October 2019, denoted as T1, and 11:00:00-18:02:39 (GPST) on 16 October 2019, denoted as T2.Applying the weighting model for the open-sky scenario, the average PDOPs amount to about 3.5 for T1, and about 3.0 for T2.The elevation-dependent and C/N 0 -dependent weighting functions, which show the best behaviours in Table 2, are used for the analysis for the two measurement scenarios.The signal standard deviations obtained with the sequential LS-VCE method (Section 2.2) and the overbounding mean and standard deviations (Section 4.1) were used for computing the positions and the protection levels, respectively.Although the rover receiver in the test of the suburban scenario is surrounded by metal infrastructure and trees as shown in the middle and right panels of Figure 1, integrity monitoring of the RTK positioning in more challenging environments, e.g., urban areas, and in kinematic mode are of interest for road transport users.These scenarios may involve more complicated and fast changing multipath environments.They are not attempted in this contribution, but are expected to be addressed in future studies.
Figure 7 shows the positioning errors in the horizontal plane and in the vertical direction for both the open-sky and the suburban scenarios.A larger error scatter can be observed for the suburban scenario.The ambiguity-fixed HPEs, computed as the square root of the squared sum of the ambiguity-fixed north and east component positioning errors, and the corresponding HPLs (Equation ( 56)) are shown in Figure 8 for T1 with higher PDOP as an example for the worse case among the two tested time periods.The ambiguity-fixed HPEs, computed as the square root of the squared sum of the north and east positioning errors, and the corresponding HPLs are shown as blue and green dots, respectively.A HAL of 0.5 m (red dots) is assumed for both scenarios.The statistics of the HPEs and the HPLs are given in Table 5 for both T1 and T2.As shown in the Table 5 and Figure 8, the FAR is achieved in the open-sky scenario for almost all time epochs for both time periods.In the suburban scenario, the FAR is achieved in about 95% of the time for T1 and almost 100% for T2.For time epochs having PAR, the majority of the ambiguities are also resolved with an average proportion of above 90% among all the ambiguities.The RMS of the HPEs are at the cm-level in both scenarios with the ambiguities fixed.Although the HPLs in the suburban scenario are higher than those in the open-sky scenario, they are generally at the sub-dm to dm-level.The ambiguity-fixed HPLs always bound the corresponding HPEs, and are always bounded by the HAL of 0.5 m.This leads to an availability of 100% for both measurement scenarios for T1 and T2.In more complicated suburban environments, larger multipath could exist in the measurements and hence bias the positioning results.For ambiguity-fixed solutions, the focus of the investigation is on the phase multipath, as the code observations are largely down-weighted, and very large code multipath is assumed to be detected by the FDE procedure.As the phase multipath should not exceed a quarter cycle [47], an extra bias of 0.25 cycle is assumed for the rover phase observations ∆φ s j (s = 1) under an elevation threshold, denoted as θs .This bias is considered in the HPLs by replacing the mDD,k in Equation (49) with mDD,k + m c,k , where m c,k is a vector containing 0.25 × λ j in range for the DD phase observations, when the rover elevation angle to satellite s (s = 1) is lower than the elevation angle threshold θs .For other observations, the entry in m c,k amounts to zero.The number of the phase observations suffering from the bias m c,k increases with the increase in θs .
Taking the time period T1 as an example, Figure 9a shows the change of the mean and maximum HPLs when changing the elevation angle threshold θs .Note that the ambiguity-fixed solutions based on the dataset used in the suburban scenario are still considered fixed.Starting from an elevation mask of 10 degrees, the mean HPL (the blue line) has increased by about 41% when increasing θs to 30 degrees.It can also be observed that the maximal HPL is still below the HAL of 0.5 m when θs is set to 30 degrees, which means that the availability is so far not influenced by the bias.Note that the change of the availability is only considered to be influenced by the change of the HPLs here, and the possible reduction of the epochs achieving FAR is not further investigated in this section.
Since the HALs for different road transport applications are not clearly set yet, in particular the very challenging autonomous driving, for applications with lower HALs, the bias m c,k could play a larger role in the achievable availabilities.In Figure 9b, the change of the availabilities is illustrated for different values of HALs selected from 0.2 m down to 0.1 m.In Figure 9b it can be observed that the availabilities with HAL of 0.2 m are similar to those with a HAL of 0.5 m, i.e., above 95%.Dramatic decrease of the availabilities occurs for applications with a higher accuracy requirement in the horizontal direction, i.e., with a HAL of 0.1 m (the blue line).In such a case, the availability is below 50% with θs higher than 21 degrees.This suggests that even for ambiguity-fixed solutions for short-baseline RTK under quiet atmosphere conditions, a clear environment with limited multipath effects might be required for applications with HAL as low as 0.1 m.

Conclusions
This contribution presents the IM for ambiguity-fixed RTK solutions under different measurement scenarios.A weighting model allowing different impacts of the elevation angles and the C/N 0 is first developed.The weighting coefficients were searched for different measurement scenarios and observation types to minimise the influences of the biases on the positional results.Next, applying these weighting models, the overbounding mean and standard deviations are determined for the weighted between-receiver observations under the open-sky and suburban scenarios based on the overbounding CDFs.
The IM is performed for RTK positioning using the ARAIM algorithm for real experiments resembling the open-sky and the suburban scenarios.After the FDE procedure, HPLs are computed for the ambiguity-fixed solutions and are compared with the HPEs and the HAL.Based on the sequential exclusion algorithm for ambiguity resolution, using GPS dual-frequency observations, FAR can be achieved in about 100% and 95% of the time for the open-sky and the suburban scenarios, respectively.Using the proposed weighting model and the mean and standard deviations derived from the overbounding EMCs, the ambiguity-fixed HPLs are at the sub-dm to dm-level.All valid HPLs are below the HAL of 0.5 m, and bound the ambiguity-fixed HPEs for both scenarios.The availabilities of the ambiguity-fixed horizontal positions reached 100% in both the open-sky and the suburban scenarios.
Investigations were also performed for the suburban scenario with larger multipath.An extra bias of a quarter cycle was added for rover phase observations below an elevation threshold.It is shown that the mean HPL could be increased by above 40% when increasing this elevation threshold from 10 to 30 degrees.For applications with high accuracy requirements, i.e., with a HAL as low as 0.1 m, a dramatic decrease of the availabilities can be observed when increasing the elevation threshold.
In more complicated measurement environments such as urban canyons, a weighting model related to both the elevation angles and the C/N 0 is also expected to behave better than that dependent solely on the elevation angles, as large multipath may appear also for signals of medium or high elevations.The bias behaviours of the DD phase and code measurements for RTK positioning in such environments, possibly not only collected from the geodetic receivers but also from low-cost receivers, will be investigated in future research.

Figure 1 .
Figure 1.Rover receivers in: (left) the open-sky scenario and (middle and right) the suburban scenario.

Figure 2 .
Figure 2. Root Mean Squares (RMS) of the double-difference (DD) residuals with respect to the elevation angles (θ 1s ) in: (a) the open-sky scenario and (b) the suburban scenario.Note that the RMS of the code residuals is divided by a factor of 100.In the bottom panels, the blue, the red and the green lines are almost overwritten by the magenta lines.

Figure 3 .
Figure 3. RMS of the DD residuals with respect to the C/N 0 (S 1s j ) in: (a) the open-sky scenario and (b) the suburban scenario.Note that the RMS of the code residuals is divided by a factor of 100.In the bottom panels, the blue and the red lines are almost overwritten by the green and the magenta lines.

Figure 4 .
Figure 4. Improvement P m in (a) phase-based and (b) code-based positional results with respect to the coefficients a φ and a p .

Figure 5 .
Figure 5. L1 and L2 residuals introducing the ground truth of the baselines and the resolved ambiguities for (a) the open-sky scenario and (b) the suburban scenario.

Figure 6 .
Figure 6.Empirical and overbounding cumulative distribution functions (CDFs) for the weighted between-receiver code residuals on L1 δp v,1 in: (a) the open-sky scenario and (b) the suburban scenario.

Figure 7 .
Figure 7. Ambiguity-fixed positioning errors in (a) the horizontal plane and (b) the vertical direction during T1 for the open-sky and the suburban scenarios.

Figure 9 .
Figure 9. Change of (a) the mean and maximal HPLs and (b) the availabilities with the elevation threshold θs for adding the phase bias term m c,k .The dataset in the suburban scenario during T1 is used for the computation.

Table 1 .
Receiver and antenna types of the two baselines used in data analysis.

Table 5 .
Statistics of the ambiguity-fixed HPEs and the HPLs for the time periods T1 and T2.The values for T1 and T2 are separated by '/'.