The local ionospheric modeling by integration ground GPS observations and satellite altimetry data

The free electrons in the ionosphere have a strong impact on the propagation of radio waves. When the signals pass through the ionosphere, both their group and phase velocity are disturbed. Several space geodetic techniques such as satellite altimetry, low Earth orbit (LEO) satellite and very long baseline interferometry (VLBI) can be used to model the total electron content. At present, the classical input data for development of ionospheric models are based on dual-frequency GPS observations, However, a major problem with this observation type is the nonuniform distribution of the terrestrial GPS reference stations with large gaps notably over the sea surface and ocean where only some single stations are located on islands, leading to lower the precision of the model over these areas. In these regions the dual-frequency satellite altimeters provide precise information about the parameters of the ionosphere. Combination of GPS and satellite altimetry observations allows making best use of the advantages of their different spatial and temporal distributions. In this study, the local ionosphere modeling was done by the combination of space geodetic observations using spherical Slepian function. The combination of the data from ground GPS observations over the western part of the USA and the altimetry mission Jason-2 was performed on the normal equation level in the least-square procedure and a least-square variance component estimation (LS-VCE) was applied to take into account the different accuracy levels of the observations. The integrated ionosphere model is more accurate and more reliable than the results derived from the ground GPS observations over the oceans.


Introduction
The Earth's ionosphere is a layer of the atmosphere at an altitude of about 60 km to 2000 km above the Earth's surface, which contains enough electrons and ions to effectively interact with the electromagnetic fields.This layer plays an important role in wireless communication due to its influence on radio propagation [Blaunstein and Christodoulou 2007].Ionospheric information such as the vertical total electron content (VTEC) can be extracted from different space geodetic techniques: global navigation satellite system, satellite altimetry missions such as TOPEX/Poseidon and Jason, very long baseline interferometry (VLBI) and low Earth orbit satellites (LEOs).Each technique has its specific characteristics influencing the derived ionosphere parameters so it can be assumed that the combined model of the ionosphere can make the best use of the advantages of each particular space geodetic technique and improve the reliability and the accuracy of the VTEC determination.The classical input data for the development of the local ionosphere map is obtained from ground-based GPS observations, however the distribution of the permanent GPS stations is not globally uniform.By combination of the GPS observations with the other data it is intended to increase the local ionosphere maps precision.
In the field of combining the different space-geodetic observations for the ionospheric modeling, several studies have been done, Todorova et al. [2007] developed the global models of the ionosphere by integration of GNSS and satellite altimetry data.Then Alizadeh et al. [2011] enhanced the model accuracy by adding radio occultation measurements (RO) from space-based GPS.For regional modeling Schmidt et al. [2008] combined GNSS, RO and satellite altimetry data and in Dettmering et al. [2011] four observation types (terrestrial GPS, space-based GPS, altimetry, and VLBI) are combined into one regional three-dimensional VTEC model.
In this study, terrestrial GPS and altimetry were integrated into one regional three-dimensional VTEC model.In this procedure, all observations were combined in one joint adjustment taking into account the systematic offsets between the observation groups.Considering the different accuracy levels of the input data in the stochastic model plays a crucial part in this procedure.In this study, the different variance factors will be estimated within a least-square variance component estimation (LS-VCE).The remainder of this paper is organized as follows.Firstly, the description of the measurements used as input data will be introduced (Section 2).Then, the general model approach and the concept for data combination will be discussed (Section 3).Finally, the VTEC model results will be presented (Section 4).

Input data
The present study is based on the data from groundbased GPS and satellite altimetry observations.Both techniques allow the observation and modelling of the ionosphere, but each of them has its specific characteristics which influence the derived ionosphere parameters.The satellite altimetry data help to compensate the insufficient GPS coverage of the oceans.The 24h observations of 30 stations belong to the International GNSS Service (IGS) network (http://sopac.ucsd.edu/)and the sampling rate of the measurements considered is 30 s.The VTEC values from the satellite altimetry Jason-2 considered in this study were obtained from the dual-frequency ionosphere delay measurements in the Geophysical Data Record (GDR).Jason-2 has taken over and continued Topex/Poseidon and Jason-1 missions since June 2008 in the frame of a cooperation among the Centre national d'études spatiales (CNES), the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT), the National Aeronautics and Space Administration (NASA) and National Oceanic and Atmospheric Administration (NOAA).It is equipped with the Poseidon-3 dual-frequency radar altimeter at 14.6 Ghz (Ku-band) and 5.3 GHz (C-band).

TEC from ground GPS
In this study a method of measuring TEC directly from the differential code delay and carrier phase measurement on both the L1 and L2 frequencies is used.For this purpose, the geometry-free linear combination of pseudo range and carrier phase measurements (also termed ionospheric observable) which is formed by subtracting simultaneous pseudorange or carrier phase observations is used [Ciraolo et al. 2007].For pseudorange and carrier phase measurements the ionospheric observable can be obtained by Equations (1) and (2) respectively: where P 1 , P 2 , U 1 and U 1 are the code and carrier phase pseudo ranges on the L1 and L2 signals with f 1 = 1575.42MHz and f 2 = 1227.60MHz frequencies, respectively.I 1 and I 2 are the ionospheric refraction delays at L1 and L2, respectively.br = c(x r p1 − x r p2 ) and Br = c(T r L1 − T r L2 ) are the code and phase inter-frequency biases (IFBs) for the receiver, bs= c(x s p1 − x s p2 ) and Bs= c(T s L1 − T s L2 ) are the code and phase differential inter-frequency biases for the satellite, f p and f L are the effects of multipath and measurement noise on the pseudorange and carrier phase, respectively, m i is the wavelength of the L i carrier phase, and N i is the integer carrier phase ambiguity.Slant total electron content (STEC) can be computed using Equations (1) or (2).Although carrier phase measurements show low noise behavior in comparison with pseudo range measurements it is not favored here for calculating STEC due to the estimation of the ambiguity term in the preprocessing.This requires the unknown ambiguities to be estimated in a pre-processing step.In order to benefit from the ambiguity-independent estimates of STECs derived from the code pseudo-ranges and the high precision of the carrier phase measurement, the pseudo-range ionospheric observables are smoothed using the "carrier to code leveling process" method [Ciraolo et al. 2007, Nohutcu et al. 2010].By combining Equations ( 1) and (2) for simultaneous observations, the following equation can be obtained: (3) Because the noise and multipath term for the carrier-phase observation is much lower than that for the pseudo-range observation, it can be neglected [Hofmann-Wellenhof et al. 2008].The interval at which no cycleslip error has occurred, leading to a constant phase ambiguity, is regarded as the continuous observational arc [Hofmann-Wellenhof et al. 2008].For cycle slip detection several testing quantities have been proposed based on various combinations of GPS observations [Seeber 2003, Hofmann-Wellenhof et al. 2008].Herein, the observation file of each station has been processed individually and a single receiver test, i.e. the combination of a phase and a code range is applied for cycle slip detection [Hofmann-Wellenhof et al. 2008].The average of the geometry-free pseudo-range and carrier phase is computed for any satellite and receiver in the continuous arc [Gao et al. 1994]: Bs br bs Br Bs br bs where n is the number of measurements in the continuous arc.By subtracting Equation (4) from Equation (2) the ambiguity term is removed: (5) P ~4 is the pseudo-range ionospheric observable smoothed by the carrier-phase ionospheric one.It will minimize the effect of multipath error so it can be neglected [Nohutcu et al. 2010].Inserting ionospheric delays into Equation ( 5), STEC can be obtained from the smoothed ionospheric observable in TEC units (TECU = 10 16 el/m 2 ) by: (6) Since in this study the ionosphere is modeled as a thin single spherical layer (Figure 1) the STEC values have to be converted into the height-independent VTEC by introducing a mapping function [Schaer 1999]: (7) where R is the mean Earth radius, z and z´are the zenith angles of the satellite at the user position and the ionospheric pierce point, and H is the mean altitude which approximately corresponds to the altitude of the maximum electron density and its height can vary between 250 and 500 km depending on latitude, season, solar and geomagnetic activity conditions [Misra and Enge 2003, Seeber 2003, Hofmann-Wellenhof et al. 2008].

TEC from satellite altimetry
The advent of the new satellite altimetry technique and the use of several different frequencies allow for estimating the ionospheric electron content.Satellite radar altimeters directly measure in nadir direction.They permanently transmit signals to Earth and receive the echo from the surface.The range between the satellite and the sea surface is determined by measuring the satellite-tosurface round-trip time of a radar pulse.Therefore, the VTEC can easily be extracted from the dual-frequency altimetry observations by converting the ionospheric correction as follows [Dumont et al. 2008 where dR is the Ku-band ionospheric range correction in meter from the Geophysical Data Record (GDR) products (iono_corr_alt_ku), and f is the Ku-band radar frequency in Hz.It is noteworthy that the TEC could then be converted to the C-band ionosphere range correction using the same formula above, and the C-band frequency of 5.3 GHz.This approximation is valid for frequencies higher than 1 GHz [Langley 1996].As the 1 Hz observations are quite noisy, the smoothing and thinning of the data is advantageous.For this purpose, a simple median filter with a length of 20 s has been applied on the altimeter ionospheric range correction [Picot et al. 2008, Dettmering et al. 2011].Figure 1 shows the sample filtered and unfiltered altimeter ionospheric range correction on the Ku-band for the Jason-2 satellite, cycle 201 and pass 137.
Because of the lower orbit altitude of the altimetry satellites, the measured VTEC are expected to be lower than the values obtained from GPS.However, several studies have shown that the TEC values obtained by the satellite altimetry are higher by about 3-4 TECU [Chelton et al. 2001, Brunini et al. 2005].Thus, in the combination stage it can be assumed that the altimetry measurements have biases, this has the effect of the plasmaspheric component (i.e. the VTEC between ~1000 and ~20,000 km height above the Earth surface [Alizadeh et al. 2011]).In this study a constant daily systematic bias for the Jason-2 satellite was estimated in the combination procedure.

Spherical Slepian functions
Spherical harmonic basis can effectively be used to represent the target function as long as the modeled area covers the whole sphere and the data is distributed regularly [Chambodut et al. 2005, Mautz et al. 2005].This can be limited by data gaps, the very high degrees with small-value coefficients, and different survey data sampling density.This function is prone to error for problems in which the data measured over the part of the Earth, as the northern hemisphere or a spherical cap, since it is no longer orthogonal over the partial area .[ Beggan et al. 2013, Sharifi andFarzaneh 2016].Therefore, the construction of a local spherical harmonic basis that is orthogonal over the studied area of the sphere is required for the study of the local modeling.In order to localize these functions into a region of interest R (target region), the optimization of a local energy criterion can be utilized.This will give a new set of functions in the sense of Slepian [1983].They are band-limited to a maximum spherical harmonic degree L and, at the same time, are spatially concentrated inside a target region.In other words, Slepian functions are a particular linear combination of the spherical harmonics.Unlike the spherical harmonics they can be sorted according to their energy concentration inside the target region [Simons et al. 2006].The spherical Slepian function can be presented as a band-limited spherical harmonic expansion: where Y lm (r) is a real spherical harmonics of degree l and order m, r is the location of a point on the surface of a unit sphere X and g lm is defined as: To maximize the spatial concentration of a bandlimited function g(r) within a region R, the ratio of the norms should be maximized as: where 0 ≤ m ≤ 1 is a measure of spatial concentration.The maximization of this concentration criterion can be achieved in the spectral domain by solving the alge-braic eigenvalue problem [Simons et al. 2006]: where the elements of (L+1) 2 × (L+1) 2 localizing kernel D: [Equation ( 13) at page bottom] are given by: and g is the (L+1) 2 dimensional vector that represents a Slepian eigenfunction expressed in the spherical harmonics, i.e: g = (g 00 ...g lm ...g LL ) T (15) This 'localization' matrix is symmetric and the subspace of maximum energy is readily obtained by eigenvalue decomposition.When the signal g(r) is local, it can be approximated using the Slepian expansion truncated at the Shannon number N [Percival and Walden 1993]: where A is the area of the region as a fraction of the full sphere.The data can be approximated, yet with very good reconstruction properties within the region by [Simons 2010]: where g(r ˆ) and d n are the spherical Slepian function and unknown coefficient, respectively.
The combination of different space-geodetic observations for the regional ionosphere modeling introduces the VTEC of two or more different space-geodetic techniques into one adjustment computation.It is obvious that the combined model is usually of higher quality because it is based on a large amount of observation data and takes advantage of the particular qualities of the input data.
The least-square adjustment (Gauss-Markov model) is applied on each set of different observations and then the individual normal equations are combined by adding the relevant N-matrices: where A i and P i are the design matrix for each observation and the weighted matrix for each observation which are uncorrelated and assumed to be known, respectively.The combination algorithm implies that the individual solutions are statistically independent populations, and that each of them provides an individual covariance matrix.Due to the various observation instruments and the functional and stochastic modeling the covariance matrices are scaled by individual factors before being introduced into the combination process.In other words, In comparison with satellite altimetry, the GNSS technique has more measurements, hence in order to increase the effect of the satellite altimetry observations in the combined model, up-weighting of the related data is needed.On the other hand, due to the higher noise of the satellite altimetry measurements, a lower weight would be implemented [Alizadeh 2013].The behavior of relative weighting is similar to scaling factor when a combination of different techniques are used.These factors (v i ) are estimated within a least-square variance component estimation (LS-VCE) [Teunissen and Amiri-Simkooei 2008].Its principal idea is to transform the Gauss-Markov model into the model of condition equations [Koch 1997 (p. 239), Bahr et al. 2007].

Results and discussion
The present study is based on the data from ground GPS and satellite altimetry observations.In order to solve STEC from ground GPS observations, the receiver IFBs were calculated using the Bernese GPS software v 5.0 and the IFB values for the satellite were obtained from the Center for Orbit Determination in Europe (CODE).The STEC and VTEC values for each observation were computed as described in Section 2. The altitude for the single-layer model was set to 400 km for the calculation of VTEC and an elevation cutoff angle of degree was used.The precise orbit files, which are provided by several IGS agencies, were in-terpolated to determine satellite positions.These TEC measurements contain ionospheric electron density information around the region above the GPS network and they are used as input data for the ionospheric modeling.Figure 2 shows the distribution of the Jason-2 observation and the GPS reference stations for August 25, 2012.
For the evaluation of the presented method, first the ionospheric map was developed by means of only the ground GPS observation.For this purpose two cases were considered, first the 2D analyses were performed, assuming that the ionosphere is static for the modeling period by neglecting the relatively small temporal variations in the ionosphere.The root-meansquare-error value for the residuals is 0.70 TECU for the reference solution.The VTEC map of the two-dimensional modeling for the mid-point of the modeling period (07:45:00) UT is plotted in Figure 3    maps shown in Figure 3 has some artifacts on the lower left corner.These artifacts may be caused by the nature of the basis functions (definition of the boundaries), associated with lack of data in this region.Second, the three-dimensional approaches based on a system of three-dimensional base functions defined as the tensor product of spherical Slepian function for longitude and latitude in an Earth-fixed reference frame, as well as polynomial B-spline functions for the time [Sharifi and Farzaneh 2013].Figure 4 shows the GPS-only VTEC map at 07:45:00 UT of day 238, 2012.The RMSE value for the residuals is 0.10 TECU for the reference solution and the optimum level of B-Spline.As can be seen, by being away from the Slepian functions domain, the results become unrealistic, especially over the ocean.The comparison between Figures 4 and 3 depicts the 3D modeling, thereby localizing the temporal variations in the ionosphere more appropriately, because in 2D modeling, it is assumed that the ionosphere has been frozen in the sun-fixed reference frame.Hence, such ionospheric modeling does not capture the small-scale and high-frequency ionosphere disturbance.Although the temporal variations have been considered in 3D modeling, but it still has low precision since there are no observations over the sea surfaces.To compensate for this deficiency, the GPS and satellite altimetry have been combined.
As mentioned before, the satellite altimetry missions provide precise information about the ionosphere over the ocean, and in the view of the fact that not many GPS stations are available there, their observations provide less accurate ionosphere maps over these regions.The combination is made by applying a least-square adjustment on each set of observations and combining the normal equations.In this procedure, a constant daily bias for the Jason-2 satellite was inserted as an additional un-known parameter in the observation equation and was estimated along with the other unknown parameters.In this study the estimated bias was −3.1 TECU.
Figure 5 illustrates the combined VTEC map at 07:45:00.It can be inferred that combining the altimetry data with the GPS observations influences the map mainly over the regions where the satellite altimetry provides information.The RMSE value for the residuals was 0.08 TECU.
To evaluate the approach and its accuracy, some of the altimetry observations were set apart from the combination procedure and were not involved in the modeling process.The number of excluded observations is selected in such a way that the accuracy of the developed model would not reduce by the confidence interval of (1−a = 95%) [Alizadeh 2013].After the development of the combined model, the VTEC values at the footprints of the raw altimetry observations, formerly eliminated from the combination, were derived from the combined model.The differences between the combined model VTEC values and the raw altimetry data were utilized for the quantification of the combined model.Hence for the altimetry satellite measurements the difference has been formed according to Alizadeh [2013]: (19) where VTEC CLIM is the combined local ionospheric model (CLIM), VTEC GDR is the raw altimetry VTEC, which is not included in CLIM, { alt , m alt are the geographic latitude and longitude of the raw altimetry observation and t is the time of the altimetry observation.Figure 6 (the red plot) presents the average of the differences in all raw observations.As stated before a daily constant offset for the altimetry measurements has been   considered within the combination procedure and is estimated along with the unknown parameters in the leastsquare procedure.This negative mean bias indicates a systematic under-estimation of VTEC delivered by Jason-2 compared to the values delivered by GPS.The blue plot in Figure 6 depicts the mean VTEC biases after removing this offset.As can be seen, the new biases vary between [−1.5, 1.5] TECU with a mean value of 0.06 TECU.This proves a good agreement of the combined model with the raw altimetry observations.

Summary and conclusion
In this study, the local ionosphere modeling has been done by the combination of space geodetic observations.The combination of data from the ground GPS observations over the western part of the USA and from the altimetry mission Jason-2 was performed on the normal equation level in the least-square procedure.The combination of the altimetry data with the ground GPS data, influenced the models mainly over the regions where satellite altimetry provided additional information.The spherical Slepian function was performed for the real set of data with the least-square variance component estimation for the weighting of the individual models.The combined model from the GPS and altimetry data had a better accuracy than the GPS-only solution over the seas.The method can be adopted for the combination with the ionospheric data from the other techniques.

Figure 1 .
Figure 1.Sample filtered and unfiltered altimeter ionospheric range correction on Ku-band for the Jason-2 satellite, cycle 201 and pass 137.

Figure 2 :
Figure 2: Distribution of Jason-2 observation and GPS reference stations for August 25, 2012.

Figure 6 .
Figure 6.Mean VTEC biases before (red plot) and after (blue plot) removing the estimated offset.