EXTENSION OF RTKLIB FOR THE CALCULATION AND VALIDATION OF PROTECTION LEVELS

: System integrity (i.e. the capability of self-monitoring) and the reliability of the positions provided need to be ensured within all safety critical applications of the GPS technology. For the sake of such applications, GPS augmentations, for example Space Based Augmentation Systems (SBAS) are to be applied to achieve the required level of integrity. SBAS provides integrity in a multi-step procedure that is laid out in the Radio Technical Commission for Aeronautics (RTCA) Minimum Operational Performance Standards (MOPS) for airborne navigation equipment using GPS. Besides integrity, SBAS also improves accuracy of positioning via broadcasting corrections to reduce the most important systematic errors on standalone positioning. To quantify integrity, the protection level is defined, which is calculated from the standard deviation of the models broadcast in SBAS. Air Navigation Service Providers, airspace users and aviation authorities need to evaluate the performance of GPS systems and their augmentations. This is a necessary step to define the conditions under which GPS systems can be operationally used and which operations can be supported. For this evaluation two proprietary software are used widely in Europe: Pegasus from Eurocontrol (Butzmühlen et al. 2001) and magicGemini from GMV. Both tools provide several functionalities such as computation of position simulating MOPS-compliant receivers and determination of GNSS augmentation attributes like accuracy, integrity, continuity and availability. RTKLIB is an open source GNSS data processing and analysis tool (Takasu, 2009). The actual version (2.4.3) of RTKLIB has SBAS augmented positioning mode, but no protection level calculation is included. There is an open source project on GitHub3, a fork of RTKLIB 2.4.2 version with an option for WAAS MOPS compliant position calculation, including protection level calculation, too. This was developed by the Houghton Associates, Inc. and tested on Cygwin platform. Their development was finished in 2014. We have merged the WAAS MOPS position calculation into the newer RTKLIB release (2.4.3 beta) and made closer integration into the original RTKLIB utility program RNX2RTKP. Our enhanced RTKLIB version is also available on GitHub4 as a fork of the original RTKLIB project of Tomoji Takasu. This enhanced version was developed and tested on Ubuntu Linux 14.04 and 16.04. Raw static and kinematic data were post-processed by our enhanced RTKLIB version. Calculated SBAS positions and protection levels were compared to the results of Pegasus and magicGemini. Although the RTCA standard defines the exact formula to calculate protection levels, the numerical results of the tested software are slightly different. Accurate tests regarding the possible sources of this kind of discrepancies were carried on in order to validate our open source solution. The


INTRODUCTION 1.1 Rationale
Due to the utmost importance attached to its security perspective, aviation is one of the most conservative technologies, being especially careful with the uptake of novel solutions.This is the reason why still today, the majority of legacy navigation equipment relies on technology dating back to the 1940s.However, industry seems to be unanimous in believing that the long-term solution within the navigation domain is that of satellite navigation, counting on the multiple constellations comprising GNSS, i.e.GPS, GLONASS and Galileo, backed up by space or ground based augmentation systems, such as WAAS, EGNOS or GAGAN.
International and European standards and guidance material promote the market uptake of these solutions.At most global level, it is the International Civil Aviation Organization, which encourages the introduction of the so-called performance based navigation solutions (PBN), in its PBN Manual (ICAO Doc 9613).Within Europe this is complemented by Regulation (EU) No 716/2014 (Pilot Common Project Regulation).This describes the PBN requirements for the 24 densest terminal manoeuvring areas in the European Union.The implementing rule is currently being drafted, this will serve as the technical basis for the legislative proposal with regard to PBN approach procedures to all instrument runway ends within Europe.
Furthermore, the ICAO Global Air Navigation Plan (ICAO Doc 9570) also treats PBN implementation as a top priority, while the SESAR (Single European Sky ATM [Air Traffic Management] Research) European ATM Master Plan puts special emphasis on the introduction of PBN as well (SESAR, 2015).This focused attention is especially needed in the short term, because it is not evident how to bridge the gap between the present situation of legacy navigation equipment and that of the long-term goal.This issue is further complicated by the 'chicken or egg' problem, or the so-called 'first mover disadvantage'.In order to reap the fruits of a fully deployed satellite based navigation network, unilaterally from airport, operator and ANSP (air navigation service provider) side, the stakeholders would ideally need to take concerted action and invest in the new technological solution unanimously.This proves difficult due to the large number of very diverse players, just as well as the highly cost sensitive environment where anyone who moves first, will inevitably encounter a temporary financial disadvantage.
Beyond this financial perspective, the new approach of satellite navigation also requires the incorporation of novel safety and security perspectives, such as the term of 'protection levels'.
Protection level is a notion introduced to be able to assign a numerical value to the reliability of navigation.The accuracy values (mean errors) resulting from the models used to reduce the effect of standard errors are calculated.Using these mean errors multiplied by protection values are protection levels defined (RTCA, 2006).If the protection level reaches the alarm limit, then satellite navigation can be seen as not fulfilling the required criteria.In theory, position error is to be below protection level at all times.If this is not the case then we talk about an integrity event (Markovits-Somogyi et al., 2017.).
The paper is structured as follows: within the next subsections of the Introduction, the ways to calculate the protection level are described and the methodology is presented.Subsequently, Section 2 comprises the results of the static measurements, while Section 3 brings along the outcomes of the kinematic measurements, obtained during a real flight trial above the regional airport of Debrecen.
It has to be noted that the data were analysed with open source software.The figures have been created with Gnuplot, while the data were uploaded into PostgreSQL database to process them, and that is where they have then been queried from.

Protection level calculation
The following equations are used to calculate the horizontal and vertical protection levels in precision approach (RTCA/DO-229C): (1) HPL=K H ⋅d major VPL=K V ⋅d U where K H =6.00  (2) where d east , d north , d EN , d U can be derived from the elements of D the variance/covariance matrix: where the i th column of the geometry matrix G T is defined as follows: (4) where El i , Az i are the elevation angle and azimuth of i th satellite.
The W −1 weight matrix inverse is a diagonal matrix with the total variances of satellites: (5) where the i th variance has four components: (6) Protection level is a notion introduced to be able to assign a numerical value to the reliability of navigation.If the protection level reaches the alarm limit defined by ICAO (International Civil Aviation Organization) standards, then satellite navigation can be seen as not fulfilling the required criteria.
Protection levels calculated from a 24 hour session of raw measurements recorded on the EGNOS monitor station operated by the Department of Geodesy and Surveying at BME (Ádám et al, 2004) are shown in Figure 1.Calculations were carried out by the upgraded version of RTKLIB with protection level calculation available on the GitHub site of Zoltán Siki (Siki, 2017).The measurements using a 1 second sampling interval in the analysed session were recorded on 28 th February 2017.
The BME EGNOS monitor station was installed in 2004 in the framework of the project EGNOS Data Collection Network (EDCN) financed by Eurocontrol.At present the station consists of a Novatel V single frequency GPS and a Glonass compatible GNSS receiver, just as well a Trimble Zephyr antenna and a PC for logging and post-processing raw measurements (Soley et al, 2004).

Software development
For the sake of carrying out the computations, an enhanced version of RTKLIB 2.4.3 beta was developed (Siki, 2017).Two sources were considered to render the protection level calculation available, the first was the official RTLIB source code on GitHub, the other was an RTKLIB 2.

Protection level comparison
The same 24 hour session of raw data was post-processed by two proprietary software: Pegasus and magicGemini.In order to compare the protection levels, the same parameters were applied, e.g. the elevation cut-off angle was set to 5 degrees.Smoothing pseudoranges by phase data was set to disabled, since this option has not yet been implemented in RTKLIB.In most of the cases only slight differences were experienced, however systematic errors seem to be omnipresent ( The minimum and maximum differences can extend to a few meters and the mean value of differences are significant.Their standard deviation remains in the range of ±0.20-0.35m.
Figure 2. Histogram of differences of the horizontal and vertical protection levels calculated by RTKLIB, Pegasus and magicGemini Although there are closed formulas to compute protection levels in RTCA standards, different values were achieved in the various computation.Two reasons may explain the differing protection levels: different software filter out particular measurements, so not the same satellites are included in the positioning solution.At the same time, it also needs to be stated that the differences of the calculated variance values are minor.

Satellites filtered out
Altogether 818,671 raw pseudorange measurements were used in positioning by RTKLIB.1217 of them were excluded by magicGemini and 1305 by Pegasus, which accounts for less than 0.16 % of the total.More or less the same measurements were filtered out by the two proprietary software.The excluded measurements have mainly be taken at lower elevation angle (Figure 3).
At the moment, there is no clear explanation available as to why these measurements were filtered out.The authors assume that a strong correlation may be found with the signal-to ratio level.In the next step, the statistical parameters of protection level differences were checked again.Only those epochs (82,803 from 86,400) were taken into account, where the very same satellites were included in the positioning.The extreme values in differences disappeared, however the trend remained the same (Table 2.).
The residual differences of the PL values may not be considered as significant, even though, as based on the standard, they are still not justifiable.The minor differences in protection levels are caused by the difference in the variances.

Variances
In formula ( 6) the four components of total variance are expressed: 1. fast and long term correction, 2. the airborne receiver error, 3. the tropospheric and 4. the ionospheric delay modelling.
Fast and long term corrections are used to decrease the effect of satellite clock offset and the errors in satellite positioning.The airborne variance contains the effect of multipath, receiver and thermal noise and the variance related to code pseudorange smoothing.For vertical ionospheric delay modelling a grid model is applied.The variance of vertical ionospheric delay is calculated as a weighted average of the variance given in the grid points.The slant delay and its variance are also calculated using a mapping function depending on the elevation angle.The variance of tropospheric delay is calculated as a constant zenith delay variance (12 cm) multiplied by the mapping function of tropospheric delay.For more details, the reader is referred to the standard (RTCA/DO-229C) or to (Ciecko and Grunwald, 2017 First the minimum and maximum values together with average and standard deviation of the differences in the fourth components of variance were checked (Table 3.).
There seem to be relatively large minimum and maximum values among the fast and long term corrections.At the same time, the majority of the differences fall between -10 and +30 cm (Table 4).In the case of fast and long term correction, the differences above 50 cm may be considered extreme.Regarding airborne variance, between RTKLIB and Pegasus calculations there are no differences, while slight differences occur between RTKLIB and magicGemini results.In 50% of the cases, the difference between RTKLIB and magicGemini values of airborne variance is -18 cm (Figure 4).It can be seen that there are not any differences in the tropospheric delay variance at all.All the three software calculated it in the perfectly same way.
However, there seem to be extremely large minimum and maximum among the difference of ionospheric delay variance (Table 3.).At the same time, the majority of the differences fall between -10 and +30 cm (Table 5).In case of the ionospheric delay, the differences above 30 cm may be considered extreme (Figure 5).The majority of the differences are observed in the north-eastern direction in case of both software.
Figure 5. Extreme large differences of ionospheric delay variance.

Flight trials in Debrecen
Having processed and analysed static data, kinematic data were collected by a GNSS receiver mounted on a general aviation aircraft.Similarly to the previous examinations, these raw data were post-processed.First, the LPV (localizer performance with vertical guidance) procedures were designed for LHDC airport (Debrecen), and then, within the framework of a flight trial these were investigated.The aim of the measurements was not only to test the flight procedures, but also to process and analyse the kinematic GNSS measurements.The flight trials were carried out on the 12th and 13th of July 2016, with the aircraft flying 6 times along the predefined flight paths.The device used for the measurements was a dual-frequency TOPCON GR3 receiver, capable of receiving both the GPS L1 and L2 and the GLONASS G1 signals.The receiver was mounted on top of the dashboard of the PA-34 220T Piper Seneca aircraft (Figure 6).The relevant standard determines explicitly and very concretely how protection levels are to be calculated.Still, the results using real life data yielded different values in case of the three different software.This difference is generally in the order of some meters.Considering the practical reliability of these numbers, it can be stated this is not significant with regard to the absolute value of the protection level.In practice, this means that the resulting values may still be used for navigational purposes, as they still meet the so-called Category I. service eligibility criteria.
However, it might be of further scientific interest to pinpoint the exact differences between these solutions, which entail the discrepancies between the different protection level calculations, sometimes even in the order of a dozen meters. 33 and 4.2 fork by the Houghton Associates, Inc. just as well on GitHub (last commit in 2014).To support our work, a new fork was created on GitHub from the official 2.4.3 beta version and the other C sources were merged into it.Some changes were made within the source code of the Houghton version.Instead of the special Cygwin version with conditional compiler directives, a more seamless integration of the two versions were carried out on Linux.The tests were run on Ubuntu 14.04 and 16.04 versions but it probably will work on any Linux distribution and Windows as well, feedback is welcome.Building RTKLIB from our source code the user shall add the -ws switch to rnx2rtp program to calculate protection levels.

Figure 3 .
Figure 3. Satellites used in positioning by RTKLIB, but not by Pegasus and magicGemini

Figure 4 .
Figure 4. Histogram of differences in airborne variance calculated by RTKLIB and magicGemini Difference in ionospheric delay correction variances

Figure 6 .
Figure 6.GNSS receiver mounted on the aircraft

Figure 8 .
Figure 8. Protection level calculated by RTKLIB during the flight trial.Measurements recorded at the permanent station in Debrecen.4. SUMMARY Vertical and horizontal protection levels are fundamental indicators of the integrity of a positioning solution.In the paper, authors have investigated positioning data from static and kinematic measurements computing them with the help of three different software, one of them open source.

Table 2 .
Statistical characteristics of protection level differences calculated by different software.Only those epochs were taken into account, where the very same satellites were included in the positioning

Table 4 .
Altogether there are 32 of such values at present.Difference in fast and long term correction variances