Tomography of the ionospheric electron density with geostatistical inversion

In relation to satellite applications like global navigation satellite systems (GNSS) and remote sensing, the electron density distribution of the ionosphere has significant influence on trans-ionospheric radio signal propagation. In this paper, we develop a novel ionospheric tomography approach providing the estimation of the electron density’s spatial covariance and based on a best linear unbiased estimator of the 3-D electron density. Therefore a non-stationary and anisotropic covariance model is set up and its parameters are determined within a maximum-likelihood approach incorporating GNSS total electron content measurements and the NeQuick model as background. As a first assessment this 3-D simple kriging approach is applied to a part of Europe. We illustrate the estimated covariance model revealing the different correlation lengths in latitude and longitude direction and its non-stationarity. Furthermore, we show promising improvements of the reconstructed electron densities compared to the background model through the validation of the ionosondes Rome, Italy (RO041), and Dourbes, Belgium (DB049), with electron density profiles for 1 day.


Introduction
The ionosphere is the upper part of the atmosphere where sufficient free electrons exist to affect the propagation of radio waves, and its morphology is mainly driven by solar radiation, particle precipitation and charge exchange.
Over the last decade, global navigation satellite system (GNSS) measurements have become one of the major tools for ionospheric sounding, enabling the derivation of the total electron content (TEC) along a satellite-to-receiver ray path.There are several activities in the ionosphere community aiming to estimate or model the ionospheric electron density based on GNSS data and other ionospheric measurements.
The International Reference Ionosphere model (IRI; see Bilitza, 2001;Bilitza and Reinisch, 2008) is an empirical model based on historical ground-and space-based data.It describes monthly averages of electron densities and temperatures in an altitude range of about 50-1500 km in the non-auroral ionosphere.Another empirical model is NeQuick (see Nava et al., 2008).It is mainly driven by the monthly average solar flux F 10.7 and ionospheric F2 peak parameters computed by the International Telecommunication Union (ITU) foF2 and M(3000)F2 models; see ITU-R (1995).However, those models represent median ionospheric behavior.Consequently the inclusion of actual ionospheric measurements is essential to update the model and hence to improve the electron density characterization.
Through the years different approaches have been developed and tested for ionospheric imaging, combining actual direct or indirect measurements with empirical or physical background models.We can identify methods modifying the coefficients of an empirical model (see Brunini et al., 2011;Galkin et al., 2012), methods updating the model towards the measurements without modification of its coefficients (see Angling and Cannon, 2004;Bust and Mitchell, 2008), meth-ods combining both (see Pezzopane et al., 2011) and approaches using physical background models and including the estimation of ionospheric drivers, such as neutral winds, in the state vector (see Schunk et al., 2004;Scherliess et al., 2009;Wang et al., 2004).
Since methods for data assimilation/ionospheric imaging were first developed, iterative methods have been used as computer resource-saving approaches to assimilate data into background models, e.g., derivatives of the algebraic reconstruction technique and the successive correction method (see Daley, 1991;Stolle et al., 2002;Gerzen and Jakowski, 2012).However, such techniques have the disadvantage that the incorporation of additional information (e.g., background and measurement error covariances), which is extremely helpful for regularization of the ill-posed inverse problem behind the ionosphere imaging, is hardly foreseen.Thus, techniques which take advantage of spatial and temporal covariance information of the ionosphere, such as optimal interpolation (OI), the 3-D and 4-D variational technique, Kalmanfilter-based approaches, and geostatistical approaches such as kriging, have been applied.In general, these methods provide a best linear unbiased estimator/predictor but differ in their mathematical frameworks and thus in their practical implementation; see e.g., Lorenc (1986).
As an example, Angling and Cannon (2004) introduced the Electron Density Assimilative Model (EDAM) that incorporates different measurements into an empirical background model by means of a Kalman filter.The majority of the input data is GPS TEC derived from the ground-based GNSS stations of the International GNSS Service (IGS).However, EDAM also deals with ionospheric radio occultation (IRO), ionosonde data and in situ electron density measurements; see Angling et al. (2008).Bust et al. (2004) developed a similar approach, the Ionospheric Data Assimilation Three-Dimensional (IDA3-D) technique based on 3-D variational data assimilation.Both EDAM and IDA3-D apply an exponential time covariance model to forecast the electron density state vector and its covariance matrix from one time step to the next.The right choice of the covariance matrix of the state vector (i.e., in this case the background covariance), the determination of the time forecast model and the appropriate choice of its parameters (for instance the correlation time) are critical to these kinds of approaches.However, until now there have only been limited publications which explicitly cover these topics.
Variograms originating from geostatistics and describing the variation between measurements depending on the distance separation are a popular tool to investigate spatial covariance; see, e.g., Wackernagel (1995); Chilès and Delfiner (2012).For the provision of vertical TEC (VTEC) and its integrity/error bounds, this method is successfully applied within the Wide Area Augmentation System (WAAS) and for the generation of global ionospheric maps (GIMs); see Blanch et al. (2003) and Pèrez (2005), respectively.In particular, both applications detrend the VTEC measurements using a background model to derive the spatial covariance of the measurements or, more specifically, the error covariance of the background.Afterwards this information is used to estimate VTEC at ionospheric grid points using ordinary kriging.
However, since electron density measurements are rarely available, especially at altitudes above the F2 layer, it is difficult to obtain the electron density's spatial covariance with variograms.In this paper, we develop an approach enabling the estimation of the electron density's spatial covariance model by means of direct and indirect ionospheric measurements.Based on this information, the electron density for arbitrary points/grids is calculated using 3-D simple kriging.

Methodology
The work flow of the approach is outlined in Fig. 1.Following the general knowledge about the ionospheric behavior, we set up a parametric spatial covariance model of the 3-D electron density.Based on the ground-based slant TEC (STEC) measurements and the NeQuick model, the unknown parameters of the spatial covariance model are derived using maximum likelihood estimation (MLE).Afterwards the electron densities of a given grid are calculated by 3-D simple kriging of linear functionals, i.e., integrals, incorporating the obtained covariance model, the NeQuick model and the STEC measurements.The subsequent sections describe each step in more detail.

Spatial covariance model of electron density
In order to establish a spatial covariance model of electron density, information about the behavior of the ionosphere is necessary.Bust et al. (2004) suggested the separation of the spatial covariance model into horizontal and vertical components to take the geometric anisotropy of the ionosphere, i.e., directionally dependent correlation lengths, into account.Yue et al. (2007) and Shim et al. (2008) confirmed this approach with the analysis of GPS and incoherent scatter radar observations revealing different correlation lengths in latitude, longitude and height direction.Furthermore, the investigations of Blanch et al. (2003) and Pèrez (2005) show that the exponential covariance model might be appropriate to describe spatial dependencies of the ionosphere.
Moreover, the non-stationarity of the ionosphere should be considered within the spatial covariance model.In other words, if we assume the electron densities N e (x i ) at arbitrary coordinates x i as a Gaussian random field (N e (x 1 ),...,N e (x n )), then the corresponding cumulative distribution functions are N (µ 1 , σ 2 1 ), . .., N (µ n , σ 2 n ) described by the expectation values µ 1 , . .., µ n and variances σ 2 1 , . .., σ 2 n , vary in space.Based on this information, we set up the following spatial covariance model of the electron density with the unknown parameter vector θ = (θ 1 , . .., θ 4 ): where x i , x j represent Earth-centered, Earth-fixed (ECEF) coordinates of the WGS84 reference ellipsoid, µ(N e (x i )), µ(N e (x j )) are the expected electron densities at the coordinates x i , x j , θ 1 is the sill parameter and c h (h h ; θ 2 , θ 3 ) and c v (h v ; θ 4 ) are the horizontal and vertical spatial covariance models, respectively.The quantities c h and c v are respectively driven by their correspondent model parameters θ 2 , . .., θ 4 and the horizontal or vertical distance, h h and h v , between two coordinates, x i and x j .The horizontal covariance model is defined as By means of the normalization, the ECEF coordinates x i and x j are projected to the unit sphere, and the influence of the height component becomes negligible.Furthermore, the assumed anisotropic correlation lengths in latitude and longitude direction are modeled by a diagonal matrix containing the parameters θ 2 and θ 3 (see Chilès and Delfiner, 2012, p. 98).Furthermore, the vertical covariance model is chosen as to be where hgt i and hgt j are the corresponding heights of x i and x j over the WGS84 reference ellipsoid.
Considering Eqs. ( 1)-( 3), it becomes clear that c h and c v describe the anisotropy of the ionosphere and the expectation values µ(N e (x i )) and µ(N e (x j )) are incorporated to take into account the ionosphere's non-stationarity.For instance, let us assume N e (x 1 ) and N e (x 2 ) to be around the ionospheric F2 peak height of 300 km, and N e (x 3 ) and N e (x 4 ) at a height of about 2000 km, where the horizontal distances h h (x 1 , x 2 ) and h h (x 3 , x 4 ) are equal.Then with Eq. ( 1) it follows that cov θ (N e (x 1 ), N e (x 2 )) is usually higher than cov θ (N e (x 3 ), N e (x 4 )).

Estimation of the spatial covariance model parameters using STEC
The spatial covariance model in Sect.2.1 depends on the parameters θ = (θ 1 , . .., θ 4 ).In order to estimate them, a background model and STEC measurements are used.Subsequently we briefly describe the calculation of ionospheric STEC measurements as well as the background model.Therefore, we particularly derive the relationship between θ and the STEC measurements and outline the MLE of θ .

Background model
As The NeQuick model is currently being developed at the International Centre for Theoretical Physics (ICTP) in Trieste, Italy, and at the University of Graz, Austria (see Hochegger et al., 2000;Radicella and Leitinger, 2001;Nava et al., 2008).It is widely used in ionospheric delay and TEC estimation for trans-ionospheric ray paths (see, e.g., Kashcheyev et al., 2012).The vertical electron density profiles of the NeQuick model are modeled by summing up five semi-Epstein layers whose shape parameters, such as peak ionization, peak height and semi-thickness, are deduced from the ITU-R (ITU Radio-Communication Sector) foF2 and M(3000)F2 models (see ITU-R, 1995).Therefore, the modeled electron density distribution inherits the spatial variances provided in the ITU-R maps via the peak ionization and peak height information.Additionally, the impact of the geomagnetic field on the ionospheric plasma density distribution is determined using a specific geomagnetic parameter called modip which is calculated from the Earth's magnetic field.The NeQuick model is driven by the solar activity level, either by the Zurich sunspot number or by the solar radio flux at 10.7 cm wave length (F 10.7 index).In the present work, we used the daily

STEC measurements
GNSS STEC measurements represent integral measurements of the electron density along a ray path s extending from a satellite position to a receiver position.By the combination of GPS dual-frequency carrier-phase (L1, L2) and code pseudorange (P1, P2) measurements, we derive the low-noise carrier-phase-relative STEC and the code-relative STEC.Subsequently, the code relative STEC is smoothed by the carrier-phase relative STEC to obtain unambiguous relative STEC measurements with a low noise level.However, the relative STEC measurements are impacted by the receiver and satellite inter-frequency biases.We use a modelassisted technique to separate the ionospheric delay (i.e., absolute STEC) and the receiver and satellite inter-frequency biases.For this, the two-dimensional Neustrelitz TEC Model (NTCM) is applied together with a mapping function based on a thin-shell ionosphere at 400 km height.For details about the absolute STEC estimation and the separation of interfrequency satellite and receiver biases we refer the reader to Jakowski et al. (2011).Based on this approach, we estimate STEC for all receiver-satellite link geometries having elevation angles equal to or greater than 10 • .

Relationship between electron density covariance and STEC measurements
We assume zero-mean Gaussian distributed and uncorrelated STEC measurement errors ε s ∼ N (0, σ 2 s ) and state the STEC measurement model as follows: where N e (s) values are the electron densities along the satellite-receiver ray path s.Since the calibration of the STEC measurements is done accordingly to Jakowski et al. (2011), the assumption of uncorrelated measurement errors is tricky.However for the purpose of the work, possible cross covariance errors are not considered.Ciraolo et al. (2007) investigated the calibration errors on experimental STEC measurements determined by GPS.He found out that the leveling of the carrier to the code measurements is mainly affected by the code multipath.Consequently, a common choice of the measurement error variance σ 2 s might be defined as dependent on the elevation angle of the satellite-to-receiver configuration assuming an increasing error budget with decreasing elevation angle.In this study, we set the minimum STEC error to 1 TECU.Considering Eqs. ( 1) and ( 4), the relationship between the spatial covariance model of the electron density and the co-variance of the STEC measurements results in , N e (r))drds + cov(ε s , ε r ) with ( 6) Assuming that the STEC measurements form a Gaussian random field − −− → STEC = (STEC s 1 , . .., STEC s n ) T with the expectation values µ = (µ(STEC s 1 ), . .., µ(STEC s n )) T and the corresponding covariance matrix ( θ ) ij := cov θ (STEC s i , STEC s j ) with i, j ∈ {1, . .., n}, the multivariate Gaussian probability density function (pdf) of the STEC measurements and is defined as where | θ | is the determinant of the covariance matrix and the expectation values µ of the STEC measurements are derived from the NeQuick model.Thus, the aim is the estimation of the parameters θ maximizing the Gaussian pdf of the STEC measurements.This maximum likelihood approach is an optimization problem and can be stated as follows: The maximization problem can be transformed into a minimization problem, for which different software package solutions exist.Within this paper, we used the Python-based software SciPy to solve the problem formulated in Eq. ( 11).
In particular, the algorithm of Powell is applied, which works iteratively and performs sequential one-dimensional minimization along each variable θ 1 , . .., θ 4 without calculating derivatives of the objective function.For more details we refer the reader to Powell (1964) and Jones et al. (2015).The initial guess for the parameter vector θ is made empirically.
We assume an electron density standard deviation of about 12 % resulting into θ 1 ≈ 0.016.Furthermore, we briefly examine the maximum horizontal distance h h between two electron densities along ray paths with ionospheric piercing points in the considered reconstruction area; see Sect.3.1.At this maximum distance, the correlation is assumed to be zero for the initial guess.Based on the investigations of Shim et al. (2008), we choose to set the initial guess for the parameters θ 2 and θ 3 to about 1.5.The parameter θ 4 controls the vertical correlation length and is set to about 300 km in agreement with the analyses of Yue et al. (2007).
Figure 2 illustrates the estimated electron density covariance models for the three different latitudes 45 • N, 50 • N and 55 • N at the Greenwich meridian at a height of 300 km.The black marked circle represents the corresponding coordinate.The correlation coefficient with its surrounding points is calculated at local times 10:00, 12:00 and 14:00, and colorcoded from blue (no correlation) to red (fully correlated).
The estimated parameters of the electron density covariance confirm the assumption of anisotropic correlation lengths in latitude and longitude.Thus, its principal behavior agrees with the TEC correlation analyses of Shim et al. (2008).Furthermore, we observe the temporal evolution of the horizontal covariance reaching its peak at 12:00 on day of year (DOY) 22 (2011) coinciding with the local time variations of the TEC correlation distances described in Yue et al. (2007).

3-D simple kriging of the electron density
Once the parameters θ of the spatial covariance model of the electron density are derived, the electron density at a WGS84 coordinate x can be estimated using simple kriging of linear functionals (i.e., integrals; see Boogaart and Drobniewski, 2002) as Consequently in order to estimate the electron density at an arbitrary WGS84 coordinate x, the product x • −1 θ ∈ R 1•n forms the weights λ = (λ 1 , . .., λ n ) T , which are used to add the difference between the GNSS-based STEC measurements and the expected STEC, in an optimal way, to the expected/modeled electron density µ[N e (x)].Moreover once the weights are calculated, the simple kriging estimation error σ 2 SK (x) at a point x is derived as (see Chilès and Delfiner, 2012, p. 153) For computational efficiency, Eq. ( 13) is extended to the dual kriging equations, enabling the estimation of N e at several WGS84 locations x 1 , . .., x m simultaneously: N e(x 1 , . .., x m ) = µ[N e (x 1 ), . .., N e (x m )] T (15) where (17) 3 Regional application

Validation scenario
In this study we apply the outlined method to a part of Europe at 40-60 • N and 30 • W-30 • E for DOY 22 (2011).We chose this region mainly for two reasons.Firstly the availability of STEC measurements is relatively good, and secondly within this region we expect better performance from the NeQuick model, which represents the background of the method, and hence an important input for the covariance and electron density estimation.DOY 22 ( 2011) is within the current maximum of solar cycle 24 but reveals quiet ionospheric conditions with a F 10.7 of 84 flux units and an average geomagnetic Kp index about 1.The STEC measurements are derived from the 1 Hz GPS L1 and L2 measurements of the International GNSS Service (IGS) ground-station network.The measurements whose corresponding ionospheric piercing points at a shell height of 400 km are within the described area are used for processing.On average, about 50 IGS stations with about 300 STEC measurements are available for a 1 s epoch; see Fig. 3. Consequently, the tomography of the ionosphere is a strongly underdetermined inverse problem with extremely limited angle geometry.Since especially the height resolution is complicated (see Garcia and Crespon, 2008) we decide to make use of STEC measurements with an elevation angle down to 10 • .
For validation, we chose the ionosondes DB049 in Dourbes, Belgium, at 50.1 • N, 4.6 • E and RO041 in Rome, Italy, at 41.9 • N, 12.5 • E; see Fig. 3.At these coordinates the height profiles of the ionospheric electron density are reconstructed and compared with the available ionosonde profiles downloaded from the Space Physics Interactive Data Resource (SPIDR).The reconstructed F2 layer characteristics in particular, in terms of NmF2 and hmF2, are validated against the measurements.For this purpose electron density profiles with a 1 km height resolution are estimated.

Preliminary results
As an example for the application of the developed ionospheric tomography, Fig. 4 illustrates the electron density layers of DOY 22 (2011), 12:00 UTC, at altitudes between 200 and 350 km.On the left-hand side, the background electron densities derived using the NeQuick model are displayed, whereas the right-hand panel depicts the electron densities calculated by the 3-D kriging.The reconstructed   Figures 5 and 6 show the electron density profiles at Dourbes, Belgium, on DOY 22 (2011) at 02:00, 12:00, 13:30 and 16:00 UTC and at Rome, Italy, on DOY 22 (2011) at 02:00, 11:45, 13:30 and 16:00 UTC.For the presented profiles an improvement of the NmF2 parameter is notable but simultaneously the peak height hmF2 is apparently not correctly reconstructed.
In order to obtain a first assessment of the 3-D kriging, we derive the F2 layer peak density and height with the 3-D kriging for DOY 22 (2011).In Fig. 7 and Table 1 the reconstructed F2 layer characteristics are validated against the model with Australian ionosonde station data.They observed that GAIM reproduces foF2 better than M(3000)F2, which is related to hmF2.For a better M(3000)F2/hmF2 reconstruction, the integration of additional ionosonde profiles and the smart handling of these data within a TEC-rich environment are noted to be crucial.
The developed 3-D kriging will be validated in more detail in future work and, in particular, the issue of the hmF2 reconstruction will be addressed.Our goal is to enhance initialization of the background model by using the ionosonde F2 layer measurements, as well as the assimilation of ionosonde electron density profiles; see, e.g., Pezzopane et al. (2013); Settimi et al. (2013).

Conclusions
The presented 3-D simple kriging of the ionospheric electron density is a novel tool for ionospheric tomography and its development is still in progress.This approach is based on the estimation of the electron density's spatial covariance, which is one of the most crucial inputs for kriging and also for different data assimilation methods.We use the relationship of this covariance to the covariance of the STEC measurements and outline the possible estimation of its parameters using the STEC measurements.Compared to the ionosonde electron density profiles, for the considered DOY 22 (2011) and locations the calculated electron density profiles show a promising gain with respect to the background model, in particular for the estimated NmF2.
In this study solely ground-based STEC measurements are incorporated.Nevertheless the approach is extendable to various ionospheric measurements such as ionosonde profiles and peak density measurements, radio occultation measurements, space-based STEC measurements and in situ measurements.In the next stage of our research we will examine this topic in order to improve the estimation of spatial covariance and electron densities.Our first effort will be the integration of the ionosonde electron density profiles, since ionosonde measurements are assumed to be the most reliable and available data type, which can provide vertical information around and below the F2 layer peak.McNamara et al. (2011) showed that in specific cases the assimilation of ionosonde data alone can yield even more accurate foF2 results than those obtained with the incorporation of GNSS TEC data in addition to the ionosonde data.
Furthermore, focus will be directed towards the inclusion of temporal information, which could be done, for instance, by developing a spatial-temporal covariance function or embedding the approach into a Kalman filter environment.Subsequently, the detailed validation will be one of the most challenging tasks.Therefore we plan a study similar to Mc-Namara et al. (2008, 2011), investigating the capability of the 3-D kriging to reconstruct the ionospheric characteristics, e.g., foF2, F2 layer thickness and M(3000)F2.Based on these results, we will refine the approach for the provision of ionospheric corrections for satellite-based radar missions in regions with dense GNSS networks.

Figure 1 .
Figure 1.Flow chart of the 3-D simple kriging approach.

Figure 3 .
Figure 3.Typical measurement geometry over a part of Europe for DOY 22 (2011) using the IGS ground-station GNSS Network: IGS ground-stations (black triangles), ionospheric piercing points of the STEC measurements (blue circles).The ionosondes RO041 and DB049 (red triangles) are used for validation.
background model an arbitrary electron density model can be considered, e.g., the NeQuick or the IRI model.
e (x i )) at a coordinate x i and the STEC along a ray path s.Additionally, the NeQuick electron density background is used within the 3-D simple kriging to stabilize the tomography of the ionospheric electron density, which presents an ill-posed and strongly underdetermined inverse problem.

Table 1 .
Comparison of the relative absolute errors of the NeQuick model and the 3-D kriging at the ionosonde station locations ofDB049  and RO041 for DOY 22 (2011), all values are given in percent.