Fault-Free Protection Level Equation for CLAS PPP-RTK and Experimental Evaluations

Centimeter level augmentation system (CLAS) of the quasi-zenith satellite system (QZSS) is the first precise point positioning-real time kinematic (PPP-RTK) augmentation system of the global navigation satellite system (GNSS), which is currently providing services for Japan. CLAS broadcasts the state-space representation of correction messages along with integrity messages regarding satellite faults and the quality index of each correction. In other GNSS augmentation systems, such as the space-based augmentation system (SBAS) of GNSS, the quality indices of correction messages are used to generate fault-free protection levels that represent a position bound containing a true user position with a probability of missed detections. Although the protection level equations are well defined for the SBAS, a protection level equation for the CLAS PPP-RTK service has not been rigorously discussed in the literature. This paper proposes a fault-free protection level equation for the PPP-RTK methods that considers the probability of correct integer ambiguity fixes in the GNSS carrier phase measurements as well as the CLAS correction quality messages. The computed protection levels with position errors were experimentally compared by processing the GNSS measurements from the GNSS Earth Observation Network (GEONET) stations in Japan and the L6 messages from the CLAS broadcast using the virtual reference station-real time kinematic (VRS-RTK) techniques. Our results, based on the GEONET dataset spanning 7 days, showed that the computed protection levels using the proposed equations were larger than the position errors for all epochs. In the dataset, the RMS errors of the CLAS VRS-RTK position were 4.6 and 14 cm in the horizontal and vertical directions, respectively, whereas the horizontal protection levels ranged from 25 cm to 2.3 m and the vertical protection levels ranged from 50 cm to 5.2 m based on fault-free integrity risk of 10−7.


Introduction
The performance of the global navigation satellite system (GNSS) is enhanced by differential GNSS techniques, including ground-and space-based augmentation systems (GBAS and SBAS) of the global positioning system (GPS) [1,2]. The GBAS and SBAS were primarily developed to guide aircraft navigation, and their design approach and operational philosophies centered on system safety. The safety levels of GBAS and SBAS are quantized as system integrity measures, and one of the important integrity measures is a protection level [3,4]. The protection level represents the bound of the true position error at the risk of a missed-detection probability. In GBAS and SBAS, a position bound is computed using the well-defined protection level equations that transform the range-domain Gaussian error distributions of each visible satellite to the position-domain Gaussian error distributions via a user-to-satellite geometry. The range-domain error distributions used in the protection level equations are determined by the GBAS and SBAS service providers and are broadcast to users as quality indicators associated with each correction message. An important characteristic of the range domain error distribution is to overbound the actual error distribution of the correction messages in order for the protection levels, computed using the quality indicators, to overbound the actual user position errors with a misseddetection probability. Consequently, studies have been conducted on various overbounding techniques based on probability density function (PDF), cumulative distribution function (CDF), and paired overbounding [5][6][7][8].
Centimeter level augmentation system (CLAS) of the quasi-zenith satellite system (QZSS) is a very recent GNSS augmentation system of Japan [9][10][11]. Unlike the GBAS and SBAS, which primarily use the GNSS code phase measurements as ranging sources, CLAS is a precise point positioning-real-time kinematic (PPP-RTK) augmentation system that allows users to resolve integer ambiguities in carrier phase measurements within several minutes by using the CLAS correction messages. Therefore, the CLAS users can use the carrier phase as a ranging source and achieve precise positioning similar to that of RTK. CLAS also provides quality indicators (also called integrity messages) for each correction message so that the users can compute the protection levels. However, the CLAS protocol currently does not specify any forms of protection-level equations, and quality message generation schemes are not well known [12,13]. Owing to lack of information, the CLAS users are not encouraged to rely on the quality messages for safe navigation.
Unfortunately, the protection level equations of GBAS or SBAS cannot be directly applied to the CLAS PPP-RTK position solutions, because GBAS or SBAS do not use the carrier phase as ranging sources. Previous studies have proposed protection-level equations for PPP and RTK that use the carrier phase as ranging sources [14][15][16][17][18][19]. These studies were based on receiver autonomous integrity monitoring (RAIM) [20], and their protection levels were derived to protect against the hard-to-detect faults in ranging measurements and cycle slips. Feng proposed a Kalman filter-based RAIM for the carrier phase [14]. In this approach, the de-correlated innovations of the Kalman filter were used to detect faults. Protection levels were computed either from the covariance matrix of the Kalman filter or from using the geometry of the satellite whose faults would be most difficult to detect. References [15,16] proposed isotropy-based protection levels (IBPL), particularly designed for the PPP position solutions. Assumptions regarding the behavior of ranging errors in terms of their size, direction, and protection levels as derived from a multivariate t-distribution of measurement errors were not made in IBPL. Ahmed et al. proposed protection levels for RTK, which were modelled in a modified form from the solution-separation RAIM methods [17,18]. An integrity risk in this approach considers the mutually exclusive events of correct and incorrect ambiguity resolutions in the least-squares ambiguity decorrelation adjustment (LAMBDA) method, which was introduced in [19]. Although these studies can be applied to derive protection levels for the CLAS PPP-RTK position solutions, their protection levels would be relatively larger because they do not use the CLAS integrity messages or the fault monitoring capability of the CLAS network.
This paper proposes a protection-level equation for the PPP-RTK services broadcasting quality messages. Because both PPP-RTK and SBAS broadcast state space representation (SSR) of the correction messages, our proposed protection level equations are based on the SBAS protection level equations and are extended to reflect virtual reference station (VRS)-RTK positioning scheme, which is a method of processing the PPP-RTK correction messages on the user side. Additionally, the integrity messages of the CLAS service are assumed to have a standard deviation of a Gaussian distribution, similar to the SBAS. This study compared the protection levels computed using the CLAS integrity messages with the VRS-RTK position errors. For the VRS-RTK process, we used CLAS L6 messages and GNSS observation measurements from the GNSS Earth Observation Network (GEONET) station in Japan. The L6 integrity messages were only used to compute the protection levels, because they were significantly inflated from the actual error distribution of the correction messages; however, the VRS-RTK process used the correction error distributions used in practice, whose standard deviation was significantly smaller than that of the integrity messages.
This paper provides an overview of the architecture of CLAS and broadcasting messages in Section 2. Section 3 discusses the proposed protection-level equations for PPP-RTK. Section 4 presents examples of broadcast integrity indices and computed horizontal and vertical protection levels with the dataset. Finally, the discussion and conclusions are presented.

Overview of the CLAS Broadcast Messages
CLAS receives the GNSS observation data from approximately 1200 GEONET stations and processes the data to generate the correction and integrity messages in a Compact SSR format. The Compact SSR messages are broadcast through the L6 band by the QZSS and are defined as RTCM 3 compatible proprietary message type 4073 [13]. The Compact SSR messages have 12 subtypes consisting of correction and integrity information, and the message types are listed in Table 1 [21]. The correction messages for the orbit, clock, code bias, and phase bias are called common corrections. The zenith tropospheric delay and slant ionospheric delay for each GNSS are referred to as local corrections. Among the message subtype IDs, the quality messages are included in subtype IDs 7, 8, 9, and 12. Subtype ID 7 provides a quality indicator for the user range accuracy (URA) of each satellite, and others provide a troposphere quality indicator as well as an SSR slant total electron content (STEC) quality indicator. The quality indicator is represented by a combination of CLASS and VALUE, the values of which range from 0 to 7, as seen in Equation (1). The SSR URA and tropospheric quality indicators were converted to a physical quantity using the following Equation [21]: The SSR STEC physical quantity was read from a table relating the SSR STEC quality indicators to the SSR STEC correction uncertainty.
The total ranging error for ith satellite can be estimated by the following equation where σ i,user is a user-specific local error, such as multipath errors. σ i,sis is a signal in space error representing the correction message uncertainty of the orbit, clock, code, and carrier biases. σ i, iono is the ionospheric delay correction uncertainty provided from the STEC correction quality indicator in total electron content unit (TECU). σ i, trop represent the troposphere delay correction uncertainty. f is the frequency of the GNSS signal and E i is the satellite elevation angle. The total ranging error is calculated in centimeters, which is why both σ i,sis and σ i, trop divided by 10 is used, and the same applies for σ i, iono . Some examples of the CLAS broadcast quality indicators are presented in Section 4.

Proposed Fault-Free Protection Level Equations for CLAS PPP-RTK Service
To develop a fault-free protection-level equation for a PPP-RTK system, the current protection-level equations of SBAS are used as the baseline equations because both SBAS and PPP-RTK broadcast the correction and quality messages in SSR formats. This section provides an overview of the SBAS fault-free protection level equation, followed by the proposed PPP-RTK protection level equation.

Fault-Free Protection Level Equations of SBAS
The SBAS L1 frequency-only protection-level equation uses the broadcast correction message uncertainties and user-defined multipath and noise uncertainties. For an individual SBAS pseudo range error, the combined range error variance for the ith satellite was constructed as follows: where σ 2 f lt,i is the variance of the fast and long-term satellite clock error corrections. σ 2 U IRE,i is the variance of the user ionosphere range error correction, and σ 2 tropo,i is the variance of the tropospheric error correction, and σ 2 air,i is the variance of the multipath error excluding multipath from ground. It should be noted that each variance is determined from a zero-mean Gaussian distribution that overbounds a non-ideal Gaussian distribution of correction errors.
The pseudo range variance is used as a weighting matrix to compute an SBAS position solution such that The direction vector from a user to satellite can be formulated as Then, the standard deviation of the position estimate uncertainty in East/North/Up (ENU) coordinates is With a missed-detection probability and its corresponding Gaussian tail value of K f f md , the vertical protection level (VPL) bounding vertical position errors is Similarly, the horizontal protection level (HPL) is computed.

Overview of Least-Squares Ambiguity Decorrelation Adjustment (LAMBDA)
The PPP-RTK service allows the use of the GNSS carrier phase for ranging measurements in a manner similar to a VRS-RTK process. A VRS-RTK process typically resolves the integer ambiguities in carrier phase measurements using the LAMBDA algorithm [22,23], and its upper bound probability of correctly fixing integer ambiguities is assessed using the integer bootstrapping method [24].
With GNSS measurements of the two receivers at a short baseline, the basic observation measurements are as follows: where ∆∇ρ and ∆∇Φ are the double-difference code and carrier-phase measurements, respectively. r is the geometric range from a user to a satellite. λ f is the f frequency wavelength of a GNSS carrier. N is integer ambiguity. ε represents the uncorrected range measurements and noise. The linearized observation equation for a set of double-difference codes and carrier-phase measurements of visible satellites is where y is a vector of the double-difference code and carrier-phase measurements. a is a vector of integer ambiguities and A is a corresponding matrix with wavelengths as elements.
B is a satellite geometry matrix, and b is the relative position vector from a VRS reference position to a GNSS receiver. The popular LAMBDA method resolves a and computes the precise b in two steps. First, an unconstrained solution of Equation (10) whereb andâ are the estimated baseline and float-integer ambiguities, respectively. Q y is the covariance matrix of ε in Equation (10). The covariance matrix ofb andâ is The second procedure of LAMBDA consists of the re-parameterization and search for a. The re-parameterization of the integer vector is implemented as follows: where Z is the decorrelation transformation matrix. Then, the re-parameterized integer ambiguity is searched with respect to the following objective function: Once an optimal integer ambiguity vector,ž, is obtained from Equation (14), the original integer ambiguity vector estimate is obtained using the inverse of the transformation such thatǎ The presumed fixed baseline vector,b iš In addition, the covariance matrix ofb is as follows: where Q ρ and Q Φ are the covariance matrices of the double-difference code and the carrierphase measurements, respectively.

Proposed Fault-Free PPP-RTK Protection Level Equation
The integer ambiguity resolution through LAMBDA is a stochastic search process, and based on the bias-free estimates of float ambiguities, its upper bound probability of correct integer fixes is [24] Prob whereâ i|I is the conditional least-squares estimate of integer ambiguity. Because the correctness of the fixed integer ambiguities may significantly affect the position errors, the probability of Equation (18) must be considered in a protection-level equation. A protection level (XPL) and given fault-free risk probability (I H0req ) can be expressed as follows: wherex and x denote the estimated position and true position, respectively, in the horizontal or vertical directions. XPL H0 refers to the horizontal or vertical protection levels. Because the position error may exceed XPL H0 with correctly fixed (CF) or incorrectly fixed (IF) integer ambiguities, Equation (19) can be expanded to two conditional probabilities, as follows [9]: where PCF is the probability of the correct fix, which can be estimated from Equation (18). PIF is the probability of incorrect fix and PCF equals to 1 − PIF.
Considering Prob x − x > XPL H0 IF = 1 as a conservative approach, the faultfree risk probability presuming that integer ambiguities are correctly fixed is Assuming that the position error follows a zero-mean Gaussian distribution, the corresponding tail value of Equation (21) was used as the K f f md factor in Equation (7).
Using the total range error of CLAS, σ i , in Equation (2), a weighting matrix for a single-differenced ranging source can be similarly constructed as in Equation (4). Because VRS-RTK uses double-difference code and carrier phase measurements, Q ρ and Q Φ in Equation (17) can be expressed as where D is the double-difference matrix.
Using the K f f md derived from Equation (21) and substituting Equation (23) into Equation (17), the fault-free horizontal and vertical protection level equations were computed as follows: (3,3) .
In Equation (23), Qb is evaluated in ENU coordinates.

Evaluation of the Proposed Protection Levels for PPP-RTK
This section discusses the GNSS observation data and parsing process of the CLAS L6 broadcast messages used in this study. Subsequently, the resultant positioning error of the VRS-RTK positioning and computed protection levels are presented.

Experimental Data Processing
To evaluate the performance of the proposed protection level, we used the GNSS observation data provided by the website of the Geospatial Information Authority of Japan [25] and the L6 data provided by the website of the QZSS [26]. CLASLIB, an opensource software tool [26], was used to extract quality messages from the L6 broadcast data and generate the GNSS measurements for the VRS-RTK process.
CLASLIB is an open-source library that parses the CLAS L6 messages, and the process is shown in Figure 1. It converts the CLAS L6 messages to the observation space representation of correction messages or provides parsed L6 message Type 4073 in comma-separated value formats. CLASLIB can also generate VRS-RTK GNSS observation data from the given GNSS navigation data and the corresponding L6 broadcast messages.
In Equation (23), is evaluated in ENU coordinates.

Evaluation of the Proposed Protection Levels for PPP-RTK
This section discusses the GNSS observation data and parsing process of the CLAS L6 broadcast messages used in this study. Subsequently, the resultant positioning error of the VRS-RTK positioning and computed protection levels are presented.

Experimental Data Processing
To evaluate the performance of the proposed protection level, we used the GNSS observation data provided by the website of the Geospatial Information Authority of Japan [25] and the L6 data provided by the website of the QZSS [26]. CLASLIB, an opensource software tool [26], was used to extract quality messages from the L6 broadcast data and generate the GNSS measurements for the VRS-RTK process.
CLASLIB is an open-source library that parses the CLAS L6 messages, and the process is shown in Figure 1. It converts the CLAS L6 messages to the observation space representation of correction messages or provides parsed L6 message Type 4073 in comma-separated value formats. CLASLIB can also generate VRS-RTK GNSS observation data from the given GNSS navigation data and the corresponding L6 broadcast messages. To evaluate the proposed protection-level equations, GNSS data were obtained from a GEONET base station located in Chiyoda City, Tokyo, Japan, as shown in Figure 2. The test site is located within the CLAS service volume. The GNSS dataset consisted of the GPS and Galileo RINEX observation and navigation files collected over 7 days from 00:00:00 on 11 April 2021, to 23:59:30 on 17 April 2021, with 30 s intervals. The base station uses a Trimble NETR9 receiver and a Trimble TPSCR.G5C antenna. The numbers of visible GPS and Galileo satellites during this period are shown in Figure 3. Because the satellite constellation exhibits a trend similar to that of a daily cycle, only the number of visible satellites in 1 day is shown. Figures 4 and 5 show the time series of the σ i,iono and σ i,tropo , respectively, extracted from the dataset. Figure 6 shows the time series of the σ i,sis for the GPS and Galileo satellites. The user receiver multipath and noise uncertainty of the To evaluate the proposed protection-level equations, GNSS data were obtained from a GEONET base station located in Chiyoda City, Tokyo, Japan, as shown in Figure 2. The test site is located within the CLAS service volume. The GNSS dataset consisted of the GPS and Galileo RINEX observation and navigation files collected over 7 days from 00:00:00 on 11 April 2021, to 23:59:30 on 17 April 2021, with 30 s intervals. The base station uses a Trimble NETR9 receiver and a Trimble TPSCR.G5C antenna. The numbers of visible GPS and Galileo satellites during this period are shown in Figure 3. Because the satellite constellation exhibits a trend similar to that of a daily cycle, only the number of visible satellites in 1 day is shown. Figures 4 and 5 show the time series of the σ i,iono and σ i,tropo , respectively, extracted from the dataset. Figure 6 shows the time series of the σ i,sis for the GPS and Galileo satellites. The user receiver multipath and noise uncertainty of the code        Figure 7 shows the resultant ENU position errors of zero-baseline VRS-RTK using dual-frequency GPS and Galileo satellites with the test dataset. To determine the fixed integers, a conventional test value ratio, above 3, was used. The RMS of the positioning errors with the fixed integers is 4, 2, and 14 cm in the East, North, and up directions, respectively. At certain epochs, our VRS-RTK software failed to resolve fixed integers due to cycle slips. For this particular test, 2.5% of the data contained a float solution.      Figures 8 and 9 show the computed protection levels and the position errors for the fixed integers cases with a fault-free risk probability (I H0req ) of 10 −7 . If the test value ratio was less than 3, that is, float solutions, no protection levels were computed. As shown in Figures 8 and 9, all of the computed HPL and VPL were larger than the horizontal positioning errors (HPE) and vertical positioning errors (VPE), respectively. The RMS of the HPE was 4.6 cm whereas that of the HPL ranged from 25 cm to 2.3 m. The RMS of the VPE was 14 cm, whereas the VPL ranged from 50 cm to 5.2 m. At certain epochs, there are large peaks in the HPL and VPL, and these peaks occur when there is a jump in quality indices. These large protection levels occurred at very small percentages such that HPL was less than 1 m in 99.4% and VPL was less than 3 m in 99.6% of the dataset, respectively.  Figure 7 shows the resultant ENU position errors of zero-baseline VRS dual-frequency GPS and Galileo satellites with the test dataset. To determin integers, a conventional test value ratio, above 3, was used. The RMS of the errors with the fixed integers is 4, 2, and 14 cm in the East, North, and up respectively. At certain epochs, our VRS-RTK software failed to resolve fixed i to cycle slips. For this particular test, 2.5% of the data contained a float solutio   Figure 7 shows the resultant ENU position errors of zero-baseline VRSdual-frequency GPS and Galileo satellites with the test dataset. To determin integers, a conventional test value ratio, above 3, was used. The RMS of the p errors with the fixed integers is 4, 2, and 14 cm in the East, North, and up respectively. At certain epochs, our VRS-RTK software failed to resolve fixed in to cycle slips. For this particular test, 2.5% of the data contained a float solution indices. These large protection levels occurred at very small percentages such that HPL was less than 1 m in 99.4% and VPL was less than 3 m in 99.6% of the dataset, respectively.

Discussion
The broadcast σ i of the correction messages in SBAS or CLAS is typically much larger than the actual correction error distribution because it is intentionally inflated to overbound the residual errors in range measurements such that the protection levels also overbound the position errors after applying the correction messages. In VRS-RTK, the uncertainty of the double-difference code and carrier phase measurement, Q ρ and Q Φ , respectively, is used to solve the float baseline and integer ambiguities. Then, the float integer ambiguities and their associated covariance matrices are used in the reparameterization and search procedure for integer ambiguities. If Q ρ and Q Φ are  indices. These large protection levels occurred at very small percentages such that HPL was less than 1 m in 99.4% and VPL was less than 3 m in 99.6% of the dataset, respectively.

Discussion
The broadcast σ i of the correction messages in SBAS or CLAS is typically much larger than the actual correction error distribution because it is intentionally inflated to overbound the residual errors in range measurements such that the protection levels also overbound the position errors after applying the correction messages. In VRS-RTK, the uncertainty of the double-difference code and carrier phase measurement, Q ρ and Q Φ , respectively, is used to solve the float baseline and integer ambiguities. Then, the float integer ambiguities and their associated covariance matrices are used in the re-

Discussion
The broadcast σ i of the correction messages in SBAS or CLAS is typically much larger than the actual correction error distribution because it is intentionally inflated to overbound the residual errors in range measurements such that the protection levels also overbound the position errors after applying the correction messages. In VRS-RTK, the uncertainty of the double-difference code and carrier phase measurement, Q ρ and Q Φ , respectively, is used to solve the float baseline and integer ambiguities. Then, the float integer ambiguities and their associated covariance matrices are used in the re-parameterization and search procedure for integer ambiguities. If Q ρ and Q Φ are constructed using the broadcasted σ i from CLAS, the LAMBDA process has a very large search space and a high chance of finding incorrect integer ambiguities. Therefore, the realistic values of Q ρ and Q Φ should be used to successfully resolve the integer ambiguities in LAMBDA.
In the test results, the RMS errors in the up direction are relatively larger than in the East and North directions. The mean value of the vertical position errors had a bias of approximately 12 cm. We also observed a similar bias from zero-baseline VRS-RTK position errors with the same dataset using an open-source software tool [28]. Therefore, it is presumed that there was a bias in the reported true antenna positions of the GNSS dataset used in the test. This issue will be investigated further.
In fact, whether the receiver moves or not, protection levels can be computed in the same way as long the PPP-RTK process is adequately implemented. However, a dynamic receiver may suffer from a frequent loss and inclusion of satellites, which may overall lower the PCF of integer ambiguities and result in incorrect integer ambiguities. The research on this issue will also be carried out as future work.

Conclusions
This paper presents a fault-free protection level equation for the CLAS PPP-RTK service that broadcasted correction quality messages. The performance of the proposed protection level equations was tested using 7 days of the GNSS observation data and the CLAS L6 messages obtained at a base station in Tokyo, Japan, which was located within the CLAS service region. The test results with 7 days of GNSS data showed that the HPL and VPL were always larger than the HPE and VPE of the zero-baseline VRS-RTK solution. Furthermore, fault-free protection-level violations were not observed. It should be noted that most HPL and VPL values were less than 1 and 3 m, respectively. Since a rail automation, which is a very demanding problem from an integrity point of view, requires Horizontal Alert Limit of 1m [29], the proposed protection level would be able to fulfill stringent integrity and availability requirements for many applications using PPP-RTK services.