Regional Ionosphere Mapping Using Zero Difference GPS Carrier Phase

Ionospheric delay, can be derived from dual frequency GNSS signals, and then converted into the Vertical Total Electron


Abstract
Ionospheric delay, can be derived from dual frequency GNSS signals, and then converted into the Vertical Total Electron Contents (VTEC) along the signal path.Various models were devised to calculate VTEC.Examples of such models are the polynomial function model and spherical harmonics analysis.A common hypothesis of these models is that they are constructed based on the assumption that the entire electron content in the ionosphere is concentrated in a single thin shell at a selected height above Earth.Due to the lack of International GNSS Service (IGS stations) over the North Africa, especially over Egypt, the current paper is presented to introduce an algorithm able to produce regional TEC maps on an hourly basis.This paper demonstrates the concept and practical examples of instantaneous ionosphere mapping, based on GPS carrier phase observations using dual frequency GPS observations over Nile Delta.The developed algorithm was based on Zero-

Introduction
Global Navigation Satellite Systems (GNSS), are considered as important costeffective tools to remotely sense the Earth's ionosphere and investigate its characteristics.This is due to the global system coverage and multiple frequency data available via a world-wide network of GPS stations.Since the ionosphere is a dispersive medium, this multiple frequency data can be used to derive the integrated measurements of electron density, which is known as TEC, along the line-of sight between a given satellite and receiver.By estimating VTEC values from a network of dual frequency GPS receivers, useful information about the ionosphere can be derived.Dual-frequency GPS receivers demonstrate the number of electrons in the ionosphere layer in a column of 1 m 2 crosssection.It is called the Slant Total Electron Content (STEC) which is extending along the ray-path of the signal between the satellite and the receiver.
Various methods were devised to calculate the ionospheric delay and STEC.These methods were based on spherical harmonic expansions in the global or regional and station scale (e.g.Schaer, 1999, andWielgosz et al., 2003a).Local methods were based on twodimensional Taylor series expansions (e.g.Komjathy, 1997, Jakobsen et al. 2010, Deng et al 2009, and Masaharu et al. 2013).A commonly used model is the grid model where the area being modeled is represented by fixed grid points in latitude and longitude.In this model, TEC values are mapped to a single ionospheric shell at fixed altitude where the maximum electron density is assumed to occur (normally a value in the range 250-400 km).
The GPS receivers use the C/A and P codes to determine the pseudo range, which is a measure of the distance between the satellite and the receiver, as the electromagnetic signal travels at the speed of light, the pseudo range can be computed by simply multiplying the time offset by the speed of light.The precision for the pseudo range measurements has been traditionally about 1% of their chip lengths, which corresponds to a precision of roughly 3 m for C/A code measurements and 0.3 m for P-code measurements (Hofmann-Wellenhof et al., 2008).GPS receiver can measure the phase of the carrier wave and track the changes in the phase but the whole number of carrier cycles that lie between the satellite and the receiver is initially unknown.To use the carrier phase as an observable for estimating ionospheric error, this unknown number of cycles or ambiguity, N, has to be determined with appropriate methods (Langley, 1998).
Ionospheric mapping is defined as a technique applying simultaneously measured TEC values to generate VTEC maps referred to a specific time epoch (Stanislawska et. al. 2000).Several studies have been performed for regional ionosphere mapping (Wielgosz et. al. 2003b;Nohutcu et. al. 2010;Salih Alcay et al. 2012).There are several groups which are capable of producing regional and/or global ionospheric maps like IGS Processing Center at the Astronomical Institute of the University of Bern, the Orbit Attitude Division of the European Space Operations Centre, and At the GPS Network and Operations Group of the Jet Propulsion Laboratory, the stations from the IGS network are using to produce regional and global ionospheric maps (Feltens et al., 1996).Regional ionospheric maps are routinely produced and made available on the Internet (Sardon et al., 1995;Jakowski et al., 1996).
Due to the lack of GPS stations over the equatorial, North Africa and Atlantic in IGS network, we suggested the proposed algorithm to be able to produce ionospheric maps over Egypt and over any other GPS station accurately.Therefore, the Ionospheric activity over Nile Delta, Egypt, was investigated using three GPS stations.The GPS data used here are from regional GPS reference observed by National Research Institute of Astronomy and Geophysics (NRIAG).
The introduced new algorithm in this study is based upon utilizing dual frequency GPS carrier phase observation because of its high precision to produce Ionospheric maps using a single GPS station (Zero Differenced) under MATLAB environment.In the present paper, Sequential Least Square Adjustment (SLSA) is considered to fix the ambiguity term to overcome the singularity in the observation equation.To test and evaluate the proposed algorithm two stations from IGS were considered, and the resulted ionosphere maps were compared with the Global Ionosphere Maps (GIMs), which is generally used as a reference.

GPS Observations
The observations of dual-frequency GPS receiver(L1 and L2) setup at any station consists of two codes and two carrier phase observations in RINEX format whichwere used for present model.The observations equations for pseudo-ranges and carrier-phase measurements can be formulated as follows (Leandro, 2009;Sedeek et al., 2017): (1) Where: P1 and P2 Pseudo-range measurements on L1 and L2 frequencies, respectively, in meter; Φ1 and Φ2 carrier-phase measurements on L1 and L2 frequencies, respectively, in meter; R the geometric distance between satellite and receiver antennas, in meters C the speed of light, in meters per second; dT and dt receiver and satellite clock errors, respectively, in seconds; T the neutral troposphere delay, in meters; I the L1 frequency ionosphere delay, in meters; Γ the factor to convert the ionospheric delay from L1 to L2 frequency, N1 and N2 carrier-phase integer ambiguities on L1 and L2 frequencies, respectively, in cycles; λ1 and λ2 carrier-phase wave length on L1 and L2 frequencies, respectively, in meters; HDr,i and HDs,i receiver and satellite pseudo-range hardware delays, respectively, in metric units, where i represents the frequency (L1 or L2); hdr,i and hds,i receiver and satellite carrier-phase hardware delays, respectively, in metric units, where i represents the frequency (L1 or L2); M1 and M2 Pseudo-range multipath on L1 and L2 frequencies, respectively, in meters; E1 and E2 Other un-modeled errors of pseudo-range measurements on L1 and L2 frequencies, respectively, in meters.pbr,i and pbs,i receiver and satellite carrier-phase initial phase bias, respectively, in metric units, where i represents the carrier frequency (L1 or L2); m1 and m2 carrier-phase multipath on L1 and L2 frequencies, respectively, in meters; e1 and e2 Other un-modeled errors of carrier-phase measurements on L1 and L2 frequencies, respectively, in meters.

Fixing Ambiguity using Sequential Least Square
Using GPS carrier-phase measurements for estimating Ionospheric delay depends upon the ability to determine the ambiguities in the number of carrier-phase cycles.Normally, the exact integer carrier-phase ambiguities are determined to estimate Ionospheric delay.
This stage gives the solution of ambiguity term.The ionospheric delay B1 which would be a check of our autonomous ionospheric delay  , that will be defined later in eq. ( 13).To fix the nature of the ambiguities, we can apply Sequential Least Square Adjustment to fix the ambiguity term to be one value of all epochs of every rise of satellite.This partitioning of the observation model is represented in following forms (Guochang Xu, 2004): Where: L is the vector of observations of the ambiguity term, A is the design matrix, X is unknown parameters vector, n is the vector of unknowns and P is weight matrix of observations with size n*n.L1 and L2 are first and second groups of observations which its number n1 and n2.The design matrix A1 and A2 and its weight matrix P1 with dimension n1*n1 and P2 with dimension n2*n2 respectively.We suppose that n1 is the first five epochs and a preliminary solution xo can be calculated as follows: xo = (A1 T P1.A1) -1 A1 T .P1.L1 (7) The change due to the additional observation set L2 is denoted as Δx: By using this system of equations, final float ambiguities values are computed.

Geometry-Free Linear Combination of GPS Observables
The geometry-free linear combination of GPS observations is classically used for ionospheric investigations and it is obtained by subtracting simultaneous pseudo range (P1-P2 or C1-P2) or carrier phase observations (Ф1-Ф2).With this combination, the satellitereceiver geometrical range and all frequency independent biases are removed (Ciraolo et al., 2007).Subtracting Eq. (3) from Eq. ( 4) the geometry-free linear combination for carrier phase observations is obtained (Sedeek et al, 2017): Where: ε(Φ1 -Φ2) is the noise term in phase equation can be neglected for the sake of simplicity, hdr  is the difference carrier-phase hardware delays bias betweenL1& L2 frequency for receiver, hds  is the difference carrier-phase hardware delays bias between L1 & L2frequency for satellite, the factor γ is the factor to convert the ionospheric delay from L1 to L2 frequency The satellite and the receiver DCBs are known to be the main error sources in the estimation of the total electron content (TEC) using GPS measurements.The TEC errors that resulted from DCBs in both the satellite and receiver can reach up to several nanoseconds (ns) (Sardon and Zarraoa 1997).So in the current study, the satellite and receiver DCB were taken from IGS products (ION files). at epoch (t) is the carrier phase Geometry free (ΦGF) after repair cycle slips.In the present paper, Melbourne-Wübbena (1985) Linear Combination is utilized to detect and repair cycle slips.The epoch-differenced in carrier phase Geometry free(ΦGF) is sufficient for estimating ionospheric delay then VTEC estimation, and therefore we use the epochdifference ∂ L4 instead of L4 itself (Araujo-Pradere et al., 2007).
Where: k is epoch number ( <k<n) and  ( ) is the ionosphere delay at epoch (t0).After applying this correction to the observation at epoch (t0 + n) at eq. ( 6) we have The epoch-carrier phase Geometry free  ( ) at eq. ( 15) includes two parts, one is due to the spatial and temporal change of the ionospheric delay from epoch (t0) to (t0 + n), and another is caused by the change of the signal paths due to the movement of the satellite ambiguity term (  −   ).

Elevation and Azimuth Angles and Ionosphere Pierce Point:
To compute elevation and azimuth angle for any satellite, the receiver position in Earth Centered Earth Fixed (ECEF) is converted to geodetic coordinate ) , , ( . Then, the satellite position coordinate ) , , ( from ECEF at the specified epoch is interpolated from the IGS final orbits.The interpolated satellite position is then transformed to a local coordinate frame, East, North, and Up (ENU) system.The transferred ENU is used to calculate elevation and azimuth angles as follows (Dahiraj, 2013 andSedeek et al, 2017): Where: E and A is elevation angle and Azimuth angle of satellite, respectively.Ionospheric Pierce Point (IPP) location can be computed by providing reference station coordinates ( , ) r r   , then the geographic latitude and longitude of IPP can be computed according to elevation and azimuth angle of satellite such as follows (Dahiraj, 2013): (18) Where:  : The offset between the IPP and the receiver; E' and E: the elevation angles at the IPP and receiver.
Where: RE: is the mean radius of the spherical Earth (6371 km) H: is the height of IPP (it is taken to be 450 km) Usually, the ionosphere assumed to be concentrated on a spherical shell of infinitesimal thickness located at altitude, for example, of 450 km above Earth's surface, i.e., forming single layer model (Rocken et al., 2000).IPP is the intersection point between the satellite receiver line-of-sight, and the ionosphere shell (Figure 1).Slant total electron content (STEC) can be translated into VTEC using Single Layer Model (SLM) (21) Where F(E) is the mapping function and E' as mentioned in equation ( 19).From equations (11 & 20) Figure ( 1): Elements of the spherical ionospheric shell model (Sedeek et al.,2017).

Evaluation of the Study
To evaluate the performance of the developed algorithm, TEC maps cover a 24 hours time period at intervals of 2 hours were estimated for observations of Day 98, 2015 for two IGS stations (ANKR, BSHM).ANKR has a geographic Longitude and Latitude (32.7583˚, 39.8875˚) and BSHM (35.02˚,32.7789˚)respectively as shown in figure (7).A geographic reference frame was used to produce the epoch-specific instantaneous regional maps of the ionosphere.Figure (2) shows the IPP distribution for ANKR station.Figures (4) and ( 5) show the regional ionosphere maps for both stations using the developed ZDPID code.We use in the current study, the satellite DCB and receiver DCB for the two stations was derived from IGS products (IONEX files).It should be noted that the IGS GIMs are computed from a Global IGS network, with observation of a combination of ~ 400 permanent GNSS stations.It is observing 4~12 satellites at 30 secretes.This means that more than 250,000 VTEC observations/hour worldwide.It uses Kalman-filter approach with shell model ionosphere with slab centered 450 km.It produces maps for every 2 hours from 180 ˚W to 180˚ E with 5˚ spatial resolution in longitude and from 87.5˚S to 87.5˚N with 2.5˚ spatial resolution in latitude (Krankowski, 2016).The IGS GIMs are provided by several analysis centers (ACs).We select a region located between 20˚~ 45˚north geographic latitude and 20˚~ 45˚ longitude.This region covers the IPPs location for most of the processed epochs.The TEC values were interpolated by Kriging method with Surfer software (V.12).We grid the studied region with space 0.5˚ × 0.5˚ to create high-resolution instantaneous TEC maps (Figure5).

Results and Discussions
The observations of Day 98 (2015) for ANKR, BSHM stations of IGS were used to evaluate the performance of the introduced ZDPID algorithm.In figures (3) and ( 4) we present the regional ionosphere maps cover a 24 hourly period at intervals of 2 hours for both stations using the developed ZDPID code with cut off angle 10˚.The TEC values have been interpolated by Kriging method using Surfer software (V.12) for the selected region located between 20˚~ 45˚N and 20˚~ 45˚ E. As it mentioned before, this region covers the IPPs location for most of the processed epochs.Figure (5) shows TEC maps estimated by GIM model of day 98 (2015) for a grid of 0.5˚ × 0.5˚ to create high-resolution instantaneous TEC maps.
As it is depicted in figure ( 5), ZDPID gives a much more detailed picture of the local ionosphere map.The IGS GIMs are a combination of TEC derivation from GPS observations, as well as different TEC modeling techniques.This also explains why the TEC derived from GIMs is very smooth over the region.In contrast to the GIMs explains the global nature of GIMs maps.But ZDPID gives real perception of the local ionosphere map.
The obtained mean TEC values every 2 hours of the two IGS stations using the ZDPID are shown in figure (6).As it is shown in the figure, the mean TEC value of IGS GIM is higher than the two stations at the same time.The difference rang between ANKR station and IGS GIM is -2.1 to 3.67 TECU and this rang is -0.29 to 3.65 TECU for BSHM station.We speculate that this may be occurred due to the global nature of IGS-GIM.IGS.Analysis centers (ACs)often use TEC representation algorithms, which result in a model resolution comparable with the whole area of the region under investigation (Schaer, 1999).Wielgosz et al., (2003b) presented an example ionosphere maps for the Ohio CORS compared to the global GIMs.The GIMs general TEC level is higher by about 3-5 TECU, as compared to the maps generated using the Kriging and Multiquadric methods.
Regional ionosphere maps have been generated for both quiet and stormy days by using Bernese 5.0 PPP model by Salih Alcay et al. (2012).They found that the biggest difference between single stations based regional model and GIMs are about 6 TEC in quiet day.Results of our model confirmed that, for regional VTEC maps are generally compatible with GIMs particularly.The GIMs results suffered from the lack of stations at some areas (e.g., over the oceans), e.g.lack of data over the equatorial, North Africa, Atlantic and in-part over equatorial and southern Pacific.This shortage of data, hamper the detection of thee equatorial anomalies (Krankowski, 2016).To overcome this shortage over Egypt and similar territories, the developed algorithm in this paper is applied to improve the temporal and spatial resolution of regional VTEC maps in Egypt.6).Mean TEC results every 2 hours of the IGS stations and IGS GIMs This will verify that TEC values could be estimated for three GPS receivers of (BORG, HELW and SAID) observed by NRIAG enclosed the Nile delta (Figure 7).BORG has a geographic Longitude and Latitude (29.5737˚, 30.86335˚),HELW (31.3434˚,29.8616˚)and SAID (32.31433˚, 31.24569˚)respectively.In the current study, the satellite DCB for the three stations was derived from IGS products (IONEX files) and the receiver DCB was calculated by (Sedeek et al., 2017).Its instantaneous regional TEC maps are drawn on a (0.5° × 0.5° grid) (Figures 8, 9, and 10).The mean values for the TEC were computed every 2 hours of the three Delta stations and IGS GIMs, the results are graphed in Figure ( 12).As it is seen in the figure, the mean TEC value of IGS GIM is higher than the three stations at the most time.The difference rang between IGS GIM and SAID station is from -1.1 to 3.69 TECU, from -1.29 to 3.27 TECU for HELW station, and from 0.2 to 4.2 TECU for BORG station.Keep in mind that the GIMs computed results of IGS do not contain any stations from North Africa, and the threshold values used in the comparison are mapped based upon interpolation techniques.
Figure (11).Mean TEC results every 2 hours of the Delta stations and IGS GIMs

Conclusions
TEC maps are needed to characterize the ionospheric behavior for satellite based positioning.Many studies have been performed for TEC mapping.In the present study a new algorithm based on Zero-differenced phase Ionospheric Delay (ZDPID) is developed.The core of this algorithm is mainly depended on computing the GPS phase ambiguity resolution model by using Sequential Least Square Adjustment.The proposed algorithm has been written using MATLAB code.The computed GIMs of IGS suffered from the lack of stations at some areas (e.g., over the oceans), e.g.lack of data over the equatorial, North Africa, Atlantic and in-part over equatorial and southern Pacific.This shortage of data, hamper the detection of the equatorial anomalies.To overcome this shortage over Egypt and similar regions, and to improve the temporal and spatial resolution of regional VTEC maps in Egypt, the new developing algorithm in this study has been applied.
To evaluate the developed algorithm, the TEC values have been estimated over two stations from the (IGS), namely ANKR and BSHM.The Global Ionosphere Maps (GIMs) were used as a threshold values for comparison with the estimated TEC of the introduced algorithm.To determine the geographical location of the regional maps, coverage circle concept has been taken into consideration.Obtained results indicated that, for regional vertical TEC maps are generally compatible with GIMs.The mean TEC value of IGS GIM found to be higher than the two stations at the same time.The difference rang between ANKR station and IGS GIM is 0.46 to 3.8 TECU and for BSHM station the difference found to be in the range from 0.1 to 3.65 TECU.ZDPID gives a much more detailed picture and real perception of the local ionosphere map from single point.On the other hand, results of the verification process are clearly shown that the developed algorithm can be successfully used for generating regional ionosphere maps.
Figure (6).Mean TEC results every 2 hours of the IGS stations and IGS GIMs

Figure ( 7
Figure (7).Location of the three Delta Stations, Egypt, and IGS stations used to estimate the regional TEC value