COMPARATIVE STUDY OF METHODS FOR CALCULATING IONOSPHERIC POINTS AND DESCRIBING THE GNSS SIGNAL PATH

: Many efforts have been done in the last decades to improve the formulation of ionospheric models based on data derived from the Global Navigation Satellite System (GNSS). Despite significant improvements in estimating the electron content of the GNSS signal path, little attention has been given to the geometric precision of ionospheric points that describe the signal path. In this work, we show a pioneer comparison about the geometric quality of the ionospheric points using distinct methods. Such analysis was carried out by calculating the GNSS signal path through three methods: a well-known geometric formulation; a new method based on linear approximations; and the used by NeQuick and recommended by the International Telecommunication Union, which was used as reference. As a result, we verified that the mean error of the well-known formulation was about 0.7 km and for the new method was at the level of 10 -11 km. Also, the proposed method has the advantage to enable the calculation of ionospheric points for GNSS signals with negative elevation angles. Once negative elevation angles derived from Radio-Occultation techniques are definitively important to improve the geometrical coverage of ionospheric modeling, the proposed technique can be useful in the development of ionospheric modeling processes.


Introduction
Ionospheric modeling processes are important tools to represent many aspects of the complex and dynamic ionosphere.In the last few decades, the Global Navigation Satellite System (GNSS) measurements have become an outstanding data source for ionospheric modeling.The ionospheric models take the dispersive propriety of the ionosphere over the GNSS signals to estimate the Total Electron Content (TEC) and describe the ionosphere in global, regional or local scales (Camargo et al., 2000;Otsuka et al., 2002;Azpilicueta et al., 2006;Mitch et al., 2013).The International GNSS Service (IGS), for example, has been using GNSS data since 1998 to provide valuable information about the global ionosphere in two-dimensional (2D) Global Ionospheric Maps (GIMs) (Schaer and Gurtner, 1998;Hernández-Pajares et al., 2009).In such maps, TEC is projected in a point located in a single layer with constant height.The projection point, called Ionospheric Pierce Point (IPP), is represented at the intersection of the satellite to receiver vector with a single layer located at 350-450 km above the Earth's surface, whereas TEC is projected on the Vertical (VTEC) of the IPPs using a specific mapping function (Schaer, 1999;Hofmann-Wellenhof et al., 2001;Rao et al., 2006).
In recent years, many efforts have been conducted for representing the ionosphere as a threedimensional (3D) layer.Multiple shell approaches (Komjathy et al., 2002), tomographic algorithms (Mitchell and Spencer, 2003;Wen et al., 2012;Prol and Camargo, 2016;Prol et al., 2017) and data assimilation procedures (Schunk et al., 2004;Hajj et al., 2004;Bust and Mitchell, 2008) are examples of ionospheric modeling algorithms that stratify the ionosphere in many layers.Such algorithms enable analysis of the vertical morphology of the ionosphere and, therefore, give a more complete specification of the ionospheric layers.However, more robust techniques are used to describe the GNSS signal path.Instead of representing TEC in a unique IPP located at 350-450 km in altitude, many ionospheric points are calculated to describe the GNSS signal path inside the 3D grid (Shukla et al., 2010;Prol and Camargo, 2015).Despite significant improvements in the last decades for VTEC interpolation on 2D maps or TEC modeling for 3D grids, little attention has been given to the geometric precision of the ionospheric points that describes the GNSS path inside the ionosphere.The ionospheric points are traditionally calculated with a simple formulation based on the elevation and azimuth angles of the GNSS signal (Skone, 2002;Matsuoka and Camargo, 2004).However, such simplification adds positional errors on the ionospheric points.In the present work, we performed experimental analysis to overview the precision of the ionospheric points using this traditional formulation.The main question that we intended to answer is "what is the magnitude of the positional error of ionospheric points calculated with the most usual method?".For that, we used as reference the recommended method by the International Telecommunication Union (ITU, 2013) to describe the GNSS signal path.
It is worth mentioning that several methods for calculating the position of ionospheric points are restricted to use GNSS signals with positive elevation angles, such as the traditional method and that one recommended by ITU.However, a significant number of GNSS signals with negative elevation angles can be derived from the presently very used Radio-Occultation (RO) technique.The RO data with negative elevation angles are definitively important to tomographic and data assimilation systems, because the incomplete geometrical coverage of the GNSS data derived from ground-based receivers make the 3D ionospheric estimation an ill-conditioned problem.Therefore, with the aim of developing a simple method for calculating the GNSS signal path for both positive and negative elevation angles, we propose and evaluate a new technique in this work.In this way, Section 2 shows the formulation of the traditional method, as well as the new technique and the method recommended by ITU.Then, in Section 3 is shown the methodology adopted for the analysis.Section 4 presents the results and discussions, and finally in Section 5 is the main conclusions of the work.

Methods to describe the GNSS signal path
Three methods are presented in this work.Section 2.1 shows the most traditional method for calculating ionospheric points, and it is being referred as Geometric method.The new method is shown in Section 2.2, which is referred as a Linear method because it is based on describing the GNSS signal using the equation of a straight line.Finally, Section 2.3 shows the most robust formulation to describe the GNSS signal path, which is recommended by ITU and used in the NeQuick ionospheric model (Nava et al., 2008).Since we used the NeQuick source code for calculating the positions according to the recommendations proposed by ITU, such method is referred here as the NeQuick method.

Geometric Method Formulation
Given the geometry of a GNSS signal tracked by a ground-based receiver station, the geographic position of an ionospheric point in the signal path is traditionally calculated by the formula presented by Matsuoka andCamargo (2004), Skone (2002) and Prol and Camargo (2015).In this method, the geographic latitude ( ip  ) and longitude ( ip  ) of the ionospheric point at a specific altitude ( ip h ) is given by a geometric relation with the azimuth and elevation angles of the GNSS signal: El are the azimuth and elevation angles of the GNSS signal in radians.Usually, the formulation presented in Equations ( 1), ( 2) and (3) are used for the calculation of IPPs at a fixed altitude of ip h = 350 km or ip h = 450 km.But it is possible to use such equations for the calculation of many ionospheric points describing the signal path.Each ionospheric point is therefore related to a distinct altitude ℎ  .Since it is a simple formulation, low computational efforts are required, which makes it an interesting formulation for many ionospheric systems that requires an intense data process.Das and Shukla (2011) and Prol and Camargo (2016), for example, use such formulation for describing the entire GNSS signal path and performing ionospheric tomographic reconstruction.However, this simplified solution requires a few geometric approximations and, therefore, intrinsic errors in the calculation of the ionospheric points are unavoidable.In Section 4 we present the magnitude of such errors.Also, the geometric method is restricted to be used for GNSS signals with positive elevation angles.In order to formulate a simple formulation that is possible to be used also for negative elevation angles, we present in the next subsection (Section 2.2) a new formulation.

Linear Method Formulation
In order to obtain a method that enables the calculation of ionospheric points for negative (and positive) elevation angles, we developed a formulation that is not dependent on the elevation angle.The alternative approach has been built in terms of cartesian coordinates by considering that the GNSS signal propagates as a straight line, which is a reasonable approximation since the ionospheric refraction includes a bending of few meters.Therefore, the coefficients that represent this straight line can be estimated.Using the linear approach, the cartesian coordinates of the ionospheric point (  4), ( 5) and ( 6): where x a , y a and z a are the angular coefficients of each straight line and r X , r Y and r Z are the cartesian coordinates of the receiver.
In the proposed new method, the fundamental interest is to obtain the angular coefficients x a , y a and z a of the straight line and the distance of the ionospheric point ip d at a specific altitude ip h from the receiver.The angular coefficients are easily obtained from Equations ( 7), ( 8) and ( 9): the distance between the receiver and the satellite.
In the case of the distance ip d between the receiver and the ionospheric point at a specific altitude ip h , an iterative procedure is adopted.A first a priori value of k ip d is adopted in the first iteration k=1.Using this initial k ip d , the cartesian coordinates of the ionospheric point are calculated using equations ( 4), ( 5) and ( 6).The cartesian coordinates for iteration k are converted into latitude ) where it is possible to obtain a small difference between of k ip ip hh  using few iterations.However, in order to validate the algorithm when the method is applied to high performance, we used one hundred iterations.
Equation ( 10) was determined by considering that the a priori distance    needs to be corrected by the difference between the expected altitude ℎ  minus the a priori altitude hh  was adopted because we consider that mostly part of the GNSS signals are closer to 45° than 90°.Therefore, less iterations are needed when we look to the whole scenario.At the end, the iterative estimated altitude 1 k ip h  will have the same value as the wished altitude ℎ  .Then, the ionospheric point position for the specific altitude ip h is obtained using 4), ( 5) and ( 6), which bring us a simple algorithm to be implemented and not dependent on the calculation of the elevation angle of the GNSS signal.

NeQuick Method Formulation
The NeQuick method was formulated based on the geometry presented in Figure 1, where P1 and P2 are the points of the receiver and satellite locations and   is the ray's perigee point.The NeQuick method is based on the calculation of the ray's perigee point, which is the point of the signal path nearest to the Earth's center.Once the   position is calculated, the entire signal path approximated as a straight line can be described by the method.with 1 r and  the radius and the zenith angle at P1.The ray perigee coordinates are calculated by Equations ( 12), ( 13) and ( 14): 11 sin( ) sin( )cos( ) cos( )sin( )cos( ) where p  is the angle between P1 and the ray perigee   :  is the azimuth of P2 seen from P1, such as in (15): and  is the Earth angle on the great circle connecting the receiver (P1) and the satellite (P2), given by Equation ( 16): cos( ) sin( )sin( ) cos( ) cos( ) cos( ) Once the ray perigee coordinates are calculated, it is possible to obtain the coordinates of any point in the signal path connecting the receiver and the satellite.To compute the geocentric coordinates of any point P (having distance  from the ray perigee   ), the following formulae ( 17) is applied: being ip h the altitude of the ionospheric point projected on the signal path.
The latitude and longitude of the ionospheric point are given by Equations ( 18) up to (22): sin( ) sin( )cos( ) cos( )sin( )cos( ) sin( ) sin( )sin( )cos( ) where the angle between the ray perigee PP and the ionospheric point is given by: tan( ) and the azimuth of a satellite as seen from ray perigee   is: 22 cos( )sin( ) sin( ) sin( ) with the following formulation to obtain the great circle angle: 2 2 2 cos( ) sin( )sin( ) cos( )cos( )cos( ) The solution of all the parameters to obtain the ray perigee and the ionospheric points based on the NeQuick method is available in the source code of NeQuick implemented in FORTRAN.In this work, such source code was used for the calculation of the ionospheric points.However, a simpler solution can be obtained online in the following website http://t-ict4d.ictp.it/nequick2/,where more details on the computation applied by NeQuick are given by ITU (2013).

Investigational Method
The Geometric, Linear and NeQuick formulations does not require information about the electron density to calculate points that describe the GNSS signal path, which means that the position of the ionospheric points is independent on the solar flux and on the ionospheric state.The main factors that affect the position of the ionospheric points are related to the geometry of the signal path.Therefore, the GNSS signal path description is affected by the geometry of the mathematical model used to describe the Earth's surface and also by the azimuth and elevation angles of the GNSS signals.In such way, the comparison between the three methods described in the previous section was performed using six ground-based receivers: ALAR (-9.75°; -36.65°, in latitude and longitude) from RBMC (Rede Brasileira de Monitoramento Contínuo) and PADO (45.41°; 11.90°), HNPT (38.58°; -76.13°),BJCO (6.38°; 2.45°), BJFS (39.61°; 115.89°) and ALIC (-23.67°;133.88°) from the International GNSS Service (IGS) network.Figure 2 shows the location of the stations, where it is possible to notice that they are located at different longitudinal sectors.
The distribution of the stations selected for this study enabled to consider distinct geometry conditions for the calculation of ionospheric points at different parts of the Earth ellipsoid.GNSS satellite positions were derived for each 15 minutes from the IGS precise ephemerides (SP3 format) recorded during the date of 25/11/2013.Then it was possible to obtain a real geometry for the elevation and azimuthal angles.For each GNSS receiver, we selected a specific GNSS satellite for the analysis with a cutoff angle of 10°.In general, we used satellites that reached a maximum daily elevation angle above 85° in order to make the analysis using satellites with a wide variety of elevation angles.For example, the GNSS satellite of PRN 20 was selected for the analysis because it provided the most varied value of elevation angles for the geometry of ALAR station.
Figure 3 shows the elevation angles for the selected geometry from ALAR, ranging from 10° up to 88°.A set of analysis was then performed with the varied conditions of elevation angles.IPPs were projected with a step size of 25 km in altitude, beginning at 50 km and ending in 2000 km.These altitudes were selected because they are commonly used to perform the integration of TEC by 3D ionospheric modeling systems, and also because 50 km and 2000 km can be considered as the bottom and top limits of the ionosphere.Figure 4 shows an example of a projected GNSS signal with an elevation angle of 10° for ALAR and PRN 20.It can be seen a bending on the GNSS signal, which is not a physical propriety.Indeed, this curvature occurs because the GNSS signal is being represented in terms of latitude, longitude and altitude, which follows the curvature of the ellipsoidal Earth.However, the signal path would be a straight line if the GNSS signal was represented in terms of cartesian coordinates.The position of the ionospheric points for each altitude was calculated using the three methods presented.The comparison between the methods was assessed by calculating the discrepancies between them and using the results from the NeQuick formulation as reference.The ionospheric points derived from NeQuick formulation were defined as reference because, as shown in the previous section, they are derived from a more robust formulation, which is also recommended by ITU for describing the GNSS signal path.The discrepancies between methods were calculated in terms of the 3D distance by using the following formulations ( 23) and ( 24): ), and the error of the linear method in terms of latitude ( , ) and longitude ( , ) in order to make clear the impact of such errors for ionospheric models, that commonly refers to the ionosphere as a layer with a specific latitude and longitude resolution.In addition, in order to outline the error magnitude and to determine which method is better on representing the ionospheric point positions, we calculated the maximum discrepancies and the Root Mean Square Error (RMSE) of both geometric and linear methods.

Results
Figure 5 shows the errors in terms of 3D distance, latitude and longitude for the same GNSS signal used in the example of Figure 4 (10° for the elevation angle).As can be noticed, the total error of the geometric method (red line) is about 0.2 km for altitudes close to 50 km and increases as the altitude elevates, until reaching 2.3 km at 2000 km.In terms of latitude and longitude, it can be verified that the maximum discrepancy is about 0.02° in latitude and 0.02° in longitude for the geometric method.On the other hand, the discrepancy obtained from the linear method (blue line) tends to zero.Similar results of Figure 5 were obtained for GNSS signals with distinct elevation angles, but Figure 5 is representative of the worst cases found in the ALAR station because it was derived using a GNSS signal with the lowest elevation angle.Considering all the signals tracked by the receiver ALAR, it can be seen in Figure 6 that the discrepancy of the geometric method tends to increase as the elevation angle gets lower.In this Figure, the discrepancy is plotted in terms of the elevation angle and for a fixed altitude, i. e., each point on the graph represents the discrepancy of the ionospheric points at a specific altitude for a signal with a specific elevation angle.We are showing the discrepancy in terms of the elevation angle for three fixed altitudes in order to present the variations of the geometric method for the commonly altitudes used to represent the ionosphere as a thin-shell (350 km and 450 km) and also for the altitude with the worst performance (2000 km).In case of the linear method, the discrepancy tends to zero for all elevation angles.As can be noticed in Figure 6, the discrepancy increases as the altitude becomes higher, showing that lower positional errors on the ionospheric points are obtained when the ionosphere is projected on a thin-layer located at 350 km than for 450 km altitude.For the geometric method, the maximum discrepancy was 0.93 km, 1.08 km and 2.33 km for the altitudes of 350 km, 450 km and 2000 km, respectively.For the linear method, the maximum discrepancy was at the level of 10 -11 km for the altitudes of 350 km, 450 km and 2000 km, showing a considerable better performance.
Despite showing details of the results of ALAR, similar behavior on the discrepancies were obtained for the other analyzed receiver stations.An overview about the discrepancy of all the receivers is presented in Figure 7, where it is shown the discrepancy for the GNSS signal related to the minimum daily elevation angle for each receiver.As it can be noticed, the discrepancy increases with altitude for the geometric method, reaching up to almost 3 km at 2000 km altitude.
Although bigger errors were found at high altitudes, there is a considerable lower concentration of electron density for points away from the density peak.Therefore, retrieving TEC with the geometric method is affected mainly by positional errors for the altitudes with larger electron densities, i. e., around 350 km and 450 km.The positional error of TEC estimation is then around 0.5 and 1 km when using the geometric method.When considering all the GNSS signals used for the stations, it can be calculated the maximum errors and the Root Mean Square Error (RMSE) for each method.Table 1 shows the resulting maximum errors and RMSE, where, in general, the total RMSE of the geometric method was around 0.7 km with a maximum of 2.38 km and, in case of the linear method, the total RMSE was 0.08 x 10 -11 km with a maximum error of 1.57 x 10 -11 km.Therefore, it is noticeable that the linear method showed a smaller discrepancy, making it more consistent with NeQuick than the geometric method.In fact, such better performance is expected due to approximations to describe the GNSS signal with geometric method, which incorporate some inconsistencies in the retrieved positions, such as the limitations pointed out by the European Aviation Safety Agency in consultation paper of the deviation request ETSO-C145c#5 that can be found in http://easa.europa.eu.The maximum error of almost 3 km presented in Table 1 for the geometric method represents a deviation of approximated 0.03° in terms of latitude and longitude.Taking into account that many ionospheric modelling systems describe the ionosphere as a grid with spatial resolutions of 0.5°, the maximum discrepancy of almost 3 km is a small value for many ionospheric modeling purposes.Therefore, it can be stated that the geometric method showed an acceptable error for ionospheric modeling purposes and are capable of being used for describing the GNSS signal path for ionospheric modeling.
It is important to mention that the error magnitude of the geometric method can be considered a problem for those monitoring small-scale irregularity structures in ionospheric studies.The irregularities that cause ionospheric scintillations in the GNSS signals have dimensions with the scale size of about 400 m and, depending on the actual ionospheric variability, the irregularities can have dimensions varying from several centimeters to a few kilometers.Therefore, errors on the ionospheric points with the magnitude of few kilometers can lead to misinterpretation of the real position of the small-scale structures.Also, such error magnitude on the ionospheric points would not be recommended for ray-tracing algorithms.The main goal of ray-tracing is to find the possible signal path in a straight line and also, the bending angle of the GNSS signal, which is caused by an ionospheric refraction from the satellite to the receiver.Since the ionospheric refraction provides a bending of few meters, an error of few kilometers due to the formulation algorithms cannot be acceptable.
In addition to show a comparative better performance than the traditional geometric method, it is possible to use the linear method for GNSS signals with negative elevation angles.Figure 8 shows an example of signal paths obtained with the linear method for negative and positive elevation angles.In this example, it can be seen that the projected ionospheric points are following the expected trajectory of the GNSS satellite, in the way that the signal paths with positive elevations are close to the receiver station and the signal path with negative elevations are far away.The elevation angle of -83° is, for example, projected almost 180° away from the receiver, showing the consistency in the representation of the ionospheric points.Therefore, since the negative elevation signals are important information to overcome the geometric limitations of the ionospheric systems using ground-based GNSS receivers, the proposed method can be considered an interesting alternative for ionospheric modeling algorithms, mainly for systems using Radio-Occultation data.

Conclusions
Two main investigations were addressed in this work.The first refers to the error magnitude on the calculation of ionospheric points using the most usual method for describing the GNSS signal path.The second inquires the possibility of developing a simple method for the calculation of ionospheric points for GNSS signals with negative elevation angles.
The error magnitude for the most used method was calculated using data collected from GNSS stations located at distinct longitudes.As a result, we verified that the geometric method presented a lower performance for points at higher altitudes, and the calculated RMSE was about 0.7 km with maximum values reaching up to 3 km at altitudes of 2000 km.In terms of ionospheric modeling, the ionospheric grids are generally obtained with a spatial resolution lower than 0.5° in latitude by 0.5° in longitude.Therefore, the error magnitude of 3 km, which is around 0.03° in the planimetric resultant, can be ignored in many applications.However, this error magnitude can be considered relevant for those application monitoring and studying small-scale structures of the ionosphere, such as the irregularities that causes ionospheric scintillations in the GNSS signals.In addition, the error magnitude shows that the traditional method is not useful for ray-tracing algorithms, where the intention is to describe the signal path in a straight line and also, the bending angle on the GNSS signal that is caused by an ionospheric refraction of few meters.
The proposed new method for calculating the GNSS signal path was evaluated with the same geometry applied to calculate the error magnitude of the most usual method.As a result, it was obtained an RMSE and a maximum error of almost zero, which is significantly better than the traditional geometric method.In addition to present a comparatively better performance, the new method enables to calculate ionospheric points for GNSS signals with negative elevation angles.The satellite signals with negative elevation angles are important for overcoming the geometric restrictions of ionospheric modeling when using GNSS Therefore, the new proposed method can be considered a good alternative for future ionospheric modeling developments using Radio-Occultation data.
receiver are given as three straight line equations, such as in Equations (

Figure 1
Figure 1 indicates the geometry involved in the computation of the coordinates of the ray perigee ( ip  ,

Figure 2 :
Figure 2: Location of the used GNSS stations.

Figure 3 :
Figure 3: Range of the elevation angles of the GNSS signals emitted from the satellite identified as PRN 20 and tracked by the GNSS receiver at ALAR.

Figure 4 :
Figure 4: Example of ionospheric points of ALAR and PRN 20 projected at every 25 km in altitude from 50 km up to 2000 km in altitude.

Figure 5 :
Figure 5: Discrepancy in terms of altitude of the geometric and linear method for a GNSS signal with a low elevation angle, station ALAR and PRN 20.

Figure 6 :
Figure 6: Discrepancy in terms of the elevation angle of the geometric and linear method for all the visible GNSS signals with a cutoff angle of 10°, station ALAR and PRN 20.

Figure 7 :
Figure 7: Discrepancy in terms of altitude of the geometric and linear method for a GNSS signal with a low elevation angle and all the analyzed stations.

Figure 8 :
Figure 8: GNSS signal paths ranging from positive to negative elevation angles calculated with the linear method for station ALAR (blue dot) and PRN 20.Each red line represents the signal path projected on the Earth and the numbers is the corresponding elevation angles.

Table 1 :
Maximum error and RMSE calculated for geometric and linear methods for all the analyzed stations.