B-spline modeling of VTEC over Turkey using GPS observations

ARTICLE IN PRESS Journal of Atmospheric and Solar-Terrestrial Physics 72 (2010) 617–624

Contents lists available at ScienceDirect

Journal of Atmospheric and Solar-Terrestrial Physics journal homepage: www.elsevier.com/locate/jastp

B-spline modeling of VTEC over Turkey using GPS observations M. Nohutcu a,1, M.O. Karslioglu a,b,n, M. Schmidt c,2 a b c

Middle East Technical University (METU), Civil Engineering Department, Geomatics Engineering Division, 06531 Ankara, Turkey METU, Department of Geodetic and Geographic Information Technologies, 06531 Ankara, Turkey Deutsches Geod¨ atisches Forschungsinstitut (DGFI), Alfons-Goppel-Str. 11, 80539 Munich, Germany

a r t i c l e in fo

abstract

Article history: Received 22 August 2009 Received in revised form 22 February 2010 Accepted 28 February 2010 Available online 4 March 2010

In this study we propose two approaches to model the Vertical Total Electron Content (VTEC) of the ionosphere with quadratic B-spline functions. For the 2-D case, VTEC is modeled in a Sun-fixed reference frame. In the 3-D approach, the 2-D model is extended to represent the temporal variations in an Earth-fixed reference frame. The localizing features of B-splines allow resolving finer structures for the regions with a sufficient number of observations by increasing the level of functions. To reduce the effects of outliers, Iteratively Re-weighted Least Squares (IRLS) with a bi-square weighting function as a robust regression algorithm is applied for parameter estimation. Another iterative method LSQR is performed for the solution of the linear systems providing a regularization effect for ill-conditioned problems. B-spline approaches are applied to real data obtained from the ground-based GPS observations over Turkey. Results are compared with the solutions of the Bernese GPS Software. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Ionosphere modeling B-splines GPS Robust regression

1. Introduction The electrons in the ionosphere affect the propagation of electromagnetic waves and this effect concerns many study areas including space-based observation systems, communication systems and space weather studies (Liu and Gao, 2003). The Global Positioning System (GPS) is also adversely affected by the ionosphere since it transmits signals in the electromagnetic spectrum. However, as it uses signals of two distinct frequencies, GPS has been widely used to monitor the electron distribution of the ionosphere for about fifteen years. Dual-frequency GPS receivers can be used to determine the number of electrons in a column of 1 m2 cross-section and extending along the ray-path of the signal between the satellite and the receiver, which is called the Slant Total Electron Content (STEC). From the ionospheric point of view, the main deficiency of the ground-based GPS measurements is their low sensitivity to the vertical structure of the atmosphere. Thus, the ionosphere is often represented by a spherical layer of infinitesimal thickness in which all the electrons are concentrated at a certain height although this approach (single layer model) introduces a modeling error due to its assumptions. Accordingly, STEC is transformed into the Vertical Total Electron Content (VTEC), which is spatially a two-dimen-

n Corresponding author at: Middle East Technical University (METU), Civil Engineering Department, Geomatics Engineering Division, 06531 Ankara, Turkey. Tel.: +90 312 2102440; fax: + 90 312 2105401. E-mail addresses: [email protected] (M. Nohutcu). [email protected] (M.O. Karslioglu), schmidt@dgfi.badw.de (M. Schmidt). 1 Tel.: + 90 312 2105472; fax: + 90 312 2105401. 2 Tel.: + 49 89 23031 1123; fax: + 49 89 23031 1240.

1364-6826/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2010.02.022

sional function depending on longitude and latitude. There are many studies to produce ionospheric maps of VTEC using GPS observations (Mannucci et al., 1998; Hernandez-Pajares et al., 1999; Schaer, 1999). An overview of ionosphere modeling based on GPS observations can be found in Brunini et al. (2004). Traditionally, VTEC is represented by spherical harmonic expansions. Spherical harmonics can be effectively used to represent the target function as long as the modeled area covers the whole sphere and the data is distributed regularly. However, the drawbacks of this method for regional applications or data of heterogeneous density have been widely discussed (Chambodut et al., 2005; Mautz et al., 2005; Schmidt et al., 2007a). VTEC can be modeled in an Earth-fixed or a Sun-fixed reference frame. The ionosphere is highly variable in an Earth-fixed reference frame due to the diurnal motion of the Earth. Thus, the models in an Earth-fixed frame should either consider the time dependency or be used instantaneously, i.e. epoch-specific. However, the ionosphere is much more stable in a Sun-fixed reference frame as the Sun is the main source for its ionization. Therefore, it can be assumed as static for a certain modeling period in a Sun-fixed reference frame (Schaer et al., 1995; Wielgosz et al., 2003). In this paper, we present two different multi-dimensional approaches, as 2-D in a Sun-fixed reference frame and 3-D in an Earth-fixed reference frame, for modeling spatio-temporal variations of VTEC. For both approaches, we propose to split VTEC into a reference part and a correction term modeled as a series expansion in terms of localizing base functions. The reference part can be obtained from a given model such as the International Reference Ionosphere (IRI) (Bilitza, 2001). However, in this study we compute it from the first level solutions of the

ARTICLE IN PRESS 618

M. Nohutcu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 72 (2010) 617–624

proposed approaches. For a regional representation with computational efficiency, we apply a multi-dimensional approach based on Euclidean quadratic B-splines and tensor-products as described, e.g., by Stollnitz et al. (1995), Lyche and Schumaker (2000), Schmidt (2007), Schmidt et al. (2007b) or Zeilhofer (2008). Since this approach is based on Euclidean theory, its application is restricted to local and regional areas. Nohutcu et al. (2008) applied the algorithm proposed by Schmidt (2007) locally by using real GPS data. The B-spline concept can be used to establish a so-called multi-resolution representation (MRR), i.e. the decomposition of the ionospheric signal into a smoother version and a number of detail signals by successive low-pass filtering. MRRs based on wavelet functions are very appropriate and efficient tools for data compression; see e.g., Schmidt (2007) or Zeilhofer (2008). The B-spline approach can be expanded in order to model the density of the ionosphere by including the height as an additional variable as described, e.g., by Schmidt (2007), Schmidt et al. (2007b), Zeilhofer (2008) or Zeilhofer et al. (2009). Note that these studies consider additional data sources, such as satellite altimetry, ionosondes or GPS receivers on Low-Earth-Orbiting (LEO) satellites, in order to overcome the low sensitivity of ground-based GPS measurements to the vertical structure of the ionosphere. The paper is organized as follows: in Section 2, the methodology for extraction of ionospheric information from GPS observations is presented. The regional representation of VTEC data with the multidimensional B-spline approach is studied in Section 3. An appropriate procedure for the estimation of the unknown series coefficients is also introduced in this section. Finally, in Section 4, the procedure is applied to real GPS data, which were collected from observation sites in Turkey.

where ^i is the carrier-phase measurement in length units for Li, li the wavelength of the Li carrier, Ni the integer carrierphase ambiguity, Br and Bs are IFBs on carrier-phase measurement. STEC can be obtained from pseudorange or carrier-phase observations by combining Eq. (1) with Eqs. (2) or (3), respectively. The noise level of carrier phase measurements is significantly lower than for pseudorange ones. However, carrier phase measurements possess ambiguity terms, which are the unknown number of whole cycles of the signal between satellite and receiver and should be estimated within a preprocessing step. Alternatively, by using the carrier-phase leveling method, the ambiguity terms can be reduced from the carrierphase ionospheric observable (Ciraolo et al., 2007) as described below. By adding Eqs. (2) and (3) for an observation, following equation can be obtained: P4 þ F4 ¼ l1 N1 l2 N2 þ Br þBs þ br þbs þ eP :

Note that the noise and multi-path term for carrier-phase observation has been neglected, as it is much lower than the one for the pseudorange observation. In Eq. (4), P4 and ^4 are available from GPS observations. The ambiguity terms N1 and N2 are constant for every continuous arc of carrier-phase observations without cycle slips. Besides, the IFB terms are stable for periods of days to months so they can be treated as constants for a continuous arc (Gao et al., 1994). Thus, Eq. (4) should provide constant or very stable results and an average value /P4 + F4Sarc can be computed for a continuous arc /P4 þ F4 Sarc ¼

A dual-frequency GPS receiver can provide code (pseudorange) and carrier phase observations for L1 and L2 signals with carrier and f2 ¼1227.60 MHz, frequencies of f1 ¼1575.42 MHz respectively. The ionospheric range delay on signals depends on the total number of electrons along the ray-path between the satellite and the receiver (STEC) and the signal frequency, and can be expressed in length units by (Liu and Gao, 2003) 40:3 fi

2

STEC:

ð1Þ

ð2Þ

where P4 is the geometry-free linear combination of pseudorange measurements, Pi the pseudorange measurement for Li, br and bs are the so-called inter-frequency biases (IFBs) on pseudorange measurement due to hardware delays of the receiver and the satellite, respectively, and ep the combination of observational noise and multi-path effects in P1 and P2. The geometry-free linear combination for carrier phase observations can be written as follows:

F4 ¼ F1 F2 ¼ I2 I1 þ l1 N1 l2 N2 þBr þ Bsþ eL ,

ð3Þ

ð5Þ

where n is the number of measurements in the continuous arc. Subtracting Eq. (3) from Eq. (5), the ambiguity terms can be eliminated P~ 4 ¼ /P4 þ F4 Sarc F4  I1 I2 þbr þ bs þ /eP Sarc eL ,

ð6Þ

where P~ 4 is the pseudorange ionospheric observable smoothed with the carrier-phase ionospheric observable. Inserting ionospheric delays from Eq. (1) into Eq. (6), STEC can be obtained in TECU (1 TECU ¼1  106 el./m2) by 2

Note, the ionospheric range delay is positive for pseudorange measurements and negative for carrier phase measurements. The geometry-free linear combination of GPS signals, which is also called the ionospheric observable, is classically used for ionospheric investigations and it is obtained by subtracting simultaneous pseudorange or carrier phase observations. With this combination, the satellite-receiver geometrical range and all frequency independent biases are removed (Ciraolo et al., 2007) P4 ¼ P1 P2 ¼ I1 I2 þ br þ bs þ ep ,

n 1X ðP4 þ F4 Þi ¼ /l1 N1 l2 N2 Sarc ni¼1

þBr þ Bsþ br þbs þ/eP Sarc ,

2. TEC determination

Ii ¼ 7

ð4Þ

STEC ¼ ðP~ 4 br bs /eP Sarc þ eL Þ

2

ðf1 f2 Þ 2

2

40:3ðf2 f1 Þ

:

ð7Þ

For two-dimensional ionosphere models, STEC values are usually converted to height independent (two-dimensional) VTEC values by introducing the so-called single layer model and corresponding mapping function. In the single layer model, all electrons in the ionosphere are assumed to be contained in a shell of infinitesimal thickness. The height of this idealized layer (H) approximately corresponds to the altitude of the maximum electron density and it is usually set to values between 350 and 450 km. The single layer model mapping function FI relates VTEC and STEC (Schaer, 1999) FI ¼

STEC 1 R with sin z0 ¼ ¼ sin z, VTEC cos z0 RþH

ð8Þ

where R is the mean Earth radius, z and z0 are the zenith angles of the satellite at the receiver and at the ionospheric pierce point (IPP), i.e. the intersection point of the receiver-to-satellite line of sight with the single layer, respectively.

ARTICLE IN PRESS M. Nohutcu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 72 (2010) 617–624

619

3. B-spline modeling of VTEC 1 k=0 value

Assume that the Earth-fixed boundary values for the rectangular area under consideration are given by [lmin, lmax] [jmin, jmax] in terms of geographic longitudes and geographic latitudes, respectively. We present two approaches for the modeling of VTEC in this area by using B-spline functions

value

k=0

ð9Þ

In 2-D modeling, the ionosphere is represented as static in a Sun-fixed reference frame for a specified time interval [tmin, tmax]. At this stage, leaving the related explanation for later, assume that the reference VTEC is computed beforehand. Thus, we have VTECVTEC ¼ DVTEC from Eq. (9). The correction term DVTEC is modeled by

DVTECðs,jÞ ¼

mJ 1 1 mJ 2 1 X X

dJ1 J2 ,k1 k2 fJ1 J2 ,k1 k2 ðs,jÞ:

ð10Þ

k1 ¼ 0 k2 ¼ 0

Here fJ1 J2 ,k1 k2 ðs,jÞ are two-dimensional scaling functions of the levels J1 and J2 with respect to the Sun-fixed longitude (s) and geographic latitude (j), mJi ¼ 2Ji +2 and dJ1 J2 ,k1 k2 are unknown scaling coefficients (Schmidt, 2007). The Sun-fixed longitude of DVTEC (or IPP) can be computed by (Dach et al., 2007) s  l þUTp,

ð11Þ

where l is the geographic longitude of IPP and UT the universal time for the observation. The 2-D scaling functions can be computed by applying the tensor product approach, i.e. separating the higher dimension scaling functions into the 1-D scaling functions

fJ1 J2 ,k1 k2 ðs,jÞ ¼ fJ1 ,k1 ðsÞ fJ2 ,k2 ðjÞ:

m Nj,k ðxÞ ¼

xtkj tkj þ m tkj

m1 Nj,k ðxÞ þ

tkj þ m þ 1 x tkj þ m þ 1 tkj þ 1

m1 Nj,k þ 1 ðxÞ:

ð13Þ

Note that when their denominators are zero, the fractions should be set to zero in Eq. (13). The recursion starts with the initial values ( ) 1 if tkj rx o tkj þ 1 0 , ð14Þ Nj,k ðxÞ ¼ 0 otherwise j where j is the level, t0j ,t1j , . . ., tmj is a sequence of non-decreasing þ2 values called knots and mj ¼2j + 2 (Schmidt, 2007). For regional modeling, endpoint-interpolation on unit interval [0, 1] is applied to avoid edge effect at the boundaries, i.e. the first

0.4 0.6 variable

0.8

1

k=1

k=2

k=3

k=4

k=5

0.5

0 0

0.2

0.4

0.6 variable

0.8

1

Fig. 1. One-dimensional B-spline scaling functions for level 1 (panel a) and level 2 (panel b). k represents shift value for each scaling function.

three knots are set to zero and the last three knots are set to one while the remaining knots are set to be equally spaced (Stollnitz et al., 1995; Lyche and Schumaker, 2000). Fig. 1 shows 1-D scaling functions for levels J¼1 and J¼2. As the figure indicates, a B-spline is compactly supported, i.e. its value is zero out of a finite range. As the level is increased the scaling functions become narrower so that finer details can be represented. However, the number of scaling coefficients (unknowns) is also increased with increased number of scaling functions. As the sequence of knots tkj implies, the variable x for the scaling function fj,k(x) takes values between 0 and 1. Hence, the coordinates s and j have to be transformed into the coordinates x and y via x¼

ssmin , smax smin

y¼

jjmin , jmax jmin

ð15Þ

where the quantities smin, smax, jmin and jmax represent the boundary coordinates of the part of the ionosphere to be modeled in the Sun-fixed system. Note that smin and smax are obtained by applying the conversion in Eq. (11) to lmin at tmin and lmax at tmax, respectively. After applying the tensor product approach (Eq. (12)) and the coordinate transformations (Eq. (15)), Eq. (10) becomes

ð12Þ

Herein fJ,k(x) is a 1-D scaling function of level J, shift k and variable x (Schmidt, 2007). In this study, the normalized quadratic B-spline N2(x) is applied as the 1-D scaling function, i.e. 2 m fj,k ðxÞ ¼ Nj,k ðxÞ. Nj,k is calculated recursively by (e.g., Stollnitz et al., 1995)

0.2

1

For both approaches, VTEC obtained in the previous section is split into two parts in this study

3.1. 2-D B-spline modeling

k=3

0.5

0

where VTEC is an approximate value or the reference and DVTEC is the correction term. VTEC can be computed from a given reference model like IRI. However, in this study we compute it from the first level B-spline solutions as described below.

k=2

0

1. 2-D modeling depending on geographic latitude and Sun-fixed longitude and 2. 3-D modeling depending on geographic latitude, geographic longitude and time.

VTEC ¼ VTEC þ DVTEC,

k=1

DVTECðx,yÞ ¼

mJ 1 1 mJ 2 1 X X

dJ1 J2 ,k1 k2 fJ1 ,k1 ðxÞ fJ2 ,k2 ðyÞ:

ð16Þ

k1 ¼ 0 k2 ¼ 0

The reference VTEC in Eq. (9) is computed with the same methodology for the correction term DVTEC computation as explained above (Eqs. (10–16)), but only taking the B-spline scaling function levels as 1 (J1 ¼J2 ¼1) and dropping the correction term, i.e. VTEC ¼ VTEC . 3.2. 3-D B-spline modeling For 3-D modeling, ionosphere is represented in an Earth-fixed reference frame depending on geographic latitude and geographic longitude, besides a time variable is included in the model as a third dimension for temporal representation. Eq. (10) is expanded to 3-D as follows:

DVTECðl,j,tÞ ¼

mJ 1 1 mJ 2 1 mJ 3 1 X X X k1 ¼ 0 k2 ¼ 0 k3 ¼ 0

dJ1 J2 J3 ,k1 k2 k3 fJ1 J2 J3 ,k1 k2 k3 ðl,j,tÞ,

ð17Þ

ARTICLE IN PRESS 620

M. Nohutcu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 72 (2010) 617–624

where fJ1 J2 J3 ,k1 k2 k3 ðl,j,tÞ are 3-D B-spline scaling functions of levels J1, J2 and J3 with respect to geographic longitude (l), geographic latitude (j) and time (t); see Schmidt et al. (2008). The 3-D scaling functions are again computed by the tensor product approach

fJ1 J2 J3 ,k1 k2 k3 ðl,j,tÞ ¼ fJ1 ,k1 ðlÞ fJ2 ,k2 ðjÞ fJ3 ,k3 ðtÞ,

ð18Þ

where the 1-D scaling functions fJ,k(x) are computed as described in Section 3.1. For 3-D modeling, the area under consideration is transformed into a cuboid and the coordinate transformations are performed as analogous to Eq. (15) including the third dimension (t) with the variable z and the boundary values tmin and tmax, i.e. the starting and ending epochs for the modeling interval. Note, that the boundary values for the longitudes are lmin and lmax for Eq. (15) as we now work in an Earth-fixed frame. Finally, the equation for the correction term is written as

DVTECðx,y,zÞ ¼

mJ 1 1 mJ 2 1 mJ 3 1 X X X

dJ1 J2 J3 ,k1 k2 k3 fJ1 ,k1 ðxÞ fJ2 ,k2 ðyÞ fJ3 ,k3 ðzÞ:

k1 ¼ 0 k2 ¼ 0 k3 ¼ 0

ð19Þ 3.3. Parameter estimation Having computed the values of the scaling functions at observation points, an observation equation is obtained for both the 2-D and the 3-D modeling (Koch, 1999) y þe ¼ Xb with DðyÞ ¼ s2 P1 ,

ð20Þ

where y is the vector of observations, e the vector of observation errors, X the coefficient matrix comprising scaling function values, b the parameter vector in terms of unknown scaling coefficients dJ1 J2 ,k1 k2 , D(y) the covariance matrix of observations, s2 an unknown variance factor and P the weight matrix of the observations. The parameter vector should include the IFB values br and bs (see Eq. (7)) if they are brought to the linear system as unknowns. However, in this study they are calculated in a preprocessing step as will be described in the next section. Note that for reference (VTEC ) computation, the observation vector y is equal to the VTEC observations, i.e. yi ¼ VTECi , and it is equal to the difference between the VTEC observations and the reference (yi ¼ VTECi VTEC i ) for the calculation of correction term (DVTEC). Least squares estimation is extensively used for the solution of a linear system of equations, e.g. as in Eq. (20). However, this method is sensitive to outliers, which may distort the estimation of the parameters if they are included in the observation data. One remedy is to detect and remove these data from the least squares estimation. Alternatively, a robust parameter estimation, which is insensitive to outliers, can be employed (Koch, 1999). In this study, a robust regression algorithm, namely Iteratively Re-weighted Least Squares (IRLS) with a bi-square weighting function is applied for parameter estimation in Eq. (20) in order to reduce the effects of outliers. As its name implies, IRLS is an iterative method where for each iteration the observations are re-weighted depending on the residuals from the previous iteration. The algorithm gives lower weights to the observations that do not fit to the model well. This structure of IRLS also reduces the effects of possible errors, e.g., due to measurement noise or single layer approach. For IRLS in this study, a bi-square weighting function is used, which is given as (Cummins and Andrews, 1995) 8" 9  2 #2 > >   > > e^ i < ^= 1 for e^ i  o ks ^ ^ , ð21Þ pB ðei Þ ¼ ks > > > > : ; 0 otherwise

where e^ i is the residual for observation i, k the tuning constant ^ an estimate for the standard which is equal to 4.685 and s deviation of the residuals which is commonly taken as s^ ¼ MAD=0:6745, where MAD is the median absolute deviation of the residuals. The robust regression algorithm starts with the computation of parameters and residuals for Eq. (20) by applying least squares estimation. For each iteration, standardized residuals are calculated as suggested by DuMouchel and O’Brien (1989) e^ i , e^ adj ¼ pffiffiffiffiffiffiffiffiffiffiffi 1hi

ð22Þ

where hi are leverages which can be computed as the diagonal elements of the hat matrix H (Koch, 1999) H ¼ XðXT XÞ1 XT ,

hi ¼ Hii :

ð23Þ

At this point, new weights for each observation are computed with the weighting function given in Eq. (21) by replacing the residuals with standardized residuals. Thus, for each iteration a new weighted linear system of equations, i.e. Eq. (20), is obtained. This system is transformed into a simplified form in which the observations are uncorrelated and have equal weights. This procedure is called homogenization in least square adjustment (Koch, 1999). The iteration steps for IRLS are repeated until the estimated parameters converge. As stated before, the values of a B-spline function are non-zero for a finite range, which defines its influence zone. A coefficient can be computed only if there are observations given in this zone. Consequently, the coefficients that are not supported with observations due to data gaps are excluded from the observation equation (Schmidt et al., 2007b). Moreover, there may exist two main difficulties for the solution of the linear system in Eq. (20). Firstly, in consequence of non-uniform data distribution, the system may become ill-conditioned. Secondly, depending on the number of observations and parameters, the size of the system may become very large. Direct regularization methods, e.g. Tikhonov’s method, are extensively used in order to stabilize the solution of ill-conditioned systems. However, these methods may become impractical due to storage problems if the system is very large, as they require great deal of computer memory (Aster et al., 2005). Schmidt (2007) uses the method of variance component estimation to determine the regularization parameter and realizes this procedure by applying a fast Monte-Carlo implementation of the iterative maximum-likelihood Variance Component Estimation; for details see Koch and Kusche (2002). Here we prefer another iterative method, namely LSQR, as a remedy for the problems that are mentioned above. LSQR is a variant of the well-known Conjugate Gradient (CG) method to be applied to the normal equations. LSQR stores only a few vectors, besides the coefficient matrix and the observation vector, and updates these vectors at each iteration. This structure of LSQR allows solving extremely large problems without matrix decomposition. Interested readers can see the algorithm for LSQR in Paige and Saunders (1982). In addition, LSQR performs a regularization effect on ill-conditioned problems where the iteration number designates the amount of regularization (Hansen, 1994). As the number of iterations is increased, the amount of regularization decreases, which indicates that the maximum regularization is achieved at the first iteration. Unfortunately, like the other iterative inverse problem solvers, LSQR does not provide the covariance matrix of model parameters. However, an approximate solution for the covariance matrix is available by modifying the LSQR algorithm as described, e.g., by Yao et al. (1999).

ARTICLE IN PRESS M. Nohutcu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 72 (2010) 617–624

A convenient method to select the appropriate value of the regularization parameter for regularization methods is the L-curve, which is an L-shaped plot of the regularized solution norm (99breg992) versus the residual norm (99Xbreg y992) in log–log scale. The corner of the L-curve can be located as an approximation to the optimal regularization parameter (Hansen, 1994). In order to select the optimal number of iterations for LSQR, the functionality of the L-curve is extended to LSQR as described by Hansen (1994), i.e. the L-curve is plotted for a discrete number of iterations. The iteration number that lies at the closest point to the corner of the L-curve is selected as the optimum iteration number. In the case that the L-curve does not represent a clear L-shape in the solution of the inverse problem, the iteration number is selected by seeking a balance between the regularized solution norm, the residual norm and the smoothness of the VTEC map.

4. Application For the application, observations of 27 ground-based dual-frequency GPS receivers located over Turkey were used. The data belong to date 26.09.2007 and have a 30 s. sampling interval. Fig. 2 demonstrates the distribution of the GPS sites that were used in the study. Observations for the stations that are marked with circles in Fig. 2 were obtained within a project that was supported by the Scientific and Technological Research Council of Turkey ¨ _ITAK). Observations for the other stations (marked with (TUB ¨ _ITAK-MAM (Marmara Research triangles) were provided by TUB Center) from their permanent networks. In order to solve STEC from observations, IFB values br and bs should be either determined within a preprocessing step or included in the parameter estimation as unknowns. We prefer to calculate receiver IFB’s with Bernese GPS Software v5.0 and obtain satellite IFB’s from CODE (Center for Orbit Determination in Europe) through the internet. The altitude for the single layer model was set to 400 km for the calculation of VTEC while a 151 elevation cut-off angle was employed for the observations. Precise orbit files comprising the cartesian coordinates of the GPS satellites at each 15 min are provided by several IGS (International GNSS Service) agencies through the internet. The satellite positions for each observation were interpolated from these precise orbits using Lagrange’s formula for polynomial interpolation in order to determine STEC directions (azimuth and zenith angles) and VTEC positions (latitude and longitude of IPP). The data corresponding to the 15:00:00–16:00:00 (UT) time slice was selected as test case since it was one of the most active periods of the ionosphere over Turkey during the selected day. We chose the rectangular region between 251 and 461 in geographic

Latitude (deg)

43

40

35 25

30

35 Longitude (deg)

40

45 46

Fig. 2. Geometry of 27 GPS stations that were used in the study.

621

longitudes and 351 and 431 in geographic latitudes as the modeling area by considering the distribution of the GPS sites. VTEC was modeled for the mentioned period in 2-D and 3-D as described in Section 3. The 2-D analyses were performed with the assumption that the ionosphere is static (or frozen) for the modeling period by neglecting the relatively small temporal variations in the ionosphere in the Sun-fixed reference frame. Therefore, the 2-D modeling results reflect the average state of the ionosphere within the related period. The 2-D solution for the levels J1 ¼1 and J2 ¼1 (B-spline levels for the Sun-fixed longitude and geographic latitude, respectively) has 16 coefficients ( ¼(21 +2)2) and this solution is also the reference (VTEC ) for the higher level solutions. The root-mean-square error value for the residuals (RMSE) is 1.1560 TECU for the reference solution. RMSE values for the DVTEC solutions for levels J1 ¼2, J2 ¼2 and levels J1 ¼3, J2 ¼ 3 are 1.1433 TECU and 1.1245 TECU with 36 and 100 coefficients, respectively. VTEC maps of the 2-D analyses for the mid-point of the modeling period (15:30:00 UT) are provided in Fig. 3. In Fig. 3, the final VTEC maps for levels J1 ¼J2 ¼2 (panel c) and J1 ¼J2 ¼3 (panel e) were obtained by adding the corresponding DVTEC solutions (panels b and d, respectively) to the reference (VTEC ) solution (panel a). While the reference solution is fairly smooth, as expected, the level of details increases as the solution levels are increased. An important point for Fig. 3 to be mentioned here is that, although the region under consideration extends between 251 and 461 in geographic longitudes, the 2-D solution considers a wider area in the Sun-fixed frame due to the diurnal motion of the Earth. For example, for the test case in this study, the Sun-fixed boundary values for longitude are 701–1061, which are computed by applying the conversion in Eq. (11) to 251 (lmin) at 15:00:00 (tmin) and 461 (lmax) at 16:00:00 (tmax), respectively. Thus, it is usual to present the VTEC maps of the 2-D solutions in the Sun-fixed reference frame. However, in order to provide a comparison between 2-D and 3-D solutions, the related parts of the VTEC maps in Sun-fixed frame, corresponding to the position of the area under consideration at 15:30:00, are cropped and presented in Earth-fixed frames in Fig. 3. For the 3-D analyses, the number of coefficients for the reference level solution at levels J1 ¼1, J2 ¼1, J3 ¼ 1 (B-spline levels for the geographic longitude, geographic latitude and time, respectively) is 64 ( ¼(21 +2)3) and its RMSE value is 0.9923 TECU. RMSE values for the DVTEC solutions for levels J1 ¼2, J2 ¼ 2, J3 ¼2 and levels J1 ¼3, J2 ¼3, J3 ¼3 are 0.9598 TECU and 0.8947 TECU with 216 and 1000 coefficients, respectively. VTEC maps of the 3-D analyses for the mid-point of the modeling period (15:30:00 UT) are shown in Fig. 4. The RMSE values indicate that the 3-D solutions fit the data better than the 2-D solutions. This result is not surprising since the 3-D solutions use the time as an extra variable to represent the temporal change in the ionosphere. On the other hand, the 2-D solutions are numerically more stable as they include far less coefficients as compared to the 3-D solutions. As in the 2-D case, higher level solutions for the 3-D modeling provide more details with smaller RMSE values. In order to assess the effect of the length of the modeling period on the results of 2-D B-spline modeling in the Sun-fixed reference frame, 2-D analyses are performed for varying modeling periods. VTEC maps for modeling periods of 15, 30, 60 and 120 min are given in Fig. 5. Periods are selected to have the same mid-point with the previous analyses, i.e. 15:30:00, so the related periods are 15:22:30–15:37:30, 15:15:00–15:45:00, 15:00:00– 16:00:00 and 14:30:00–16:30:00 in UT. For a clear demonstration of the results, only the VTEC maps for levels J1 ¼2 and J2 ¼2 are given in Fig. 5, excluding the plots for references and corrections.

ARTICLE IN PRESS 622

M. Nohutcu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 72 (2010) 617–624

Latitude (deg)

15 42 10

40 38

5

36 25

30

35

0 45 TECU

40

2

15 42

42 1 40

10

40

0 38

38

5

-1 36 25

36

30

35

40

-2 45 TECU

25

30

35

40

0 45 TECU 15

2 42

42 1

10

40

40 0

38

38

5

-1 36

36 25

30

35

40

-2 45 TECU

25

30

35

40

0 45 TECU

Longitude (deg) Fig. 3. 2-D VTEC model results for 26.09.2007, 15:30 (UT). Modeling period is 15:00:00–16:00:00 (UT). (a) VTEC solution for J1 ¼ J2 ¼1 with 16 coefficients. RMSE¼ 1.1560 TECU; (b) DVTEC solution for J1 ¼J2 ¼2 with 36 coefficients. RMSE¼ 1.1433 TECU; (c) VTEC solution for J1 ¼J2 ¼ 2 by adding levels 2, 2 corrections (panel b) to reference (panel a); (d) DVTEC solution for J1 ¼J2 ¼ 3 with 100 coefficients. RMSE¼ 1.1245 TECU; (e) VTEC solution for J1 ¼ J2 ¼3 by adding levels 3, 3 corrections (panel d) to reference (panel a).

15

Latitude (deg)

42 10

40 38

5

36 25

30

35

40

0 45 TECU 15

42

2

40

42 10

40 0

38

38 -2

36 25

30

35

40

45 TECU

5

36 25

30

35

40

0 45 TECU

15 42

2

42 10

40

40 0

38

38 -2

36 25

30

35

40

45 TECU

5

36 25

30

35

40

0 45 TECU

Longitude (deg) Fig. 4. 3-D VTEC model results for 26.09.2007, 15:30 (UT). Modeling period is 15:00:00–16:00:00 (UT). (a) VTEC solution for J1 ¼J2 ¼ J3 ¼ 1 with 64 coefficients. RMSE¼ 0.9923 TECU; (b) DVTEC solution for J1 ¼ J2 ¼J3 ¼ 2 with 216 coefficients. RMSE¼ 0.9598 TECU; (c) VTEC solution for J1 ¼J2 ¼ J3 ¼2 by adding levels 2, 2, 2 corrections (panel b) to reference (panel a); (d) DVTEC solution for J1 ¼J2 ¼ J3 ¼3 with 1000 coefficients. RMSE ¼0.8947 TECU; (e) VTEC solution for J1 ¼J2 ¼J3 ¼ 3 by adding levels 3, 3, 3 corrections (panel d) to reference (panel a).

The RMSE values for the periods of 15, 30, 60 and 120 min are 1.0886, 1.1168, 1.1433 and 1.1708 TECU, respectively. The increase in RMSE values for longer modeling periods is an

expected result as the 2-D model of VTEC in a Sun-fixed reference frame demonstrates the average state of the ionosphere for the related period as mentioned before. Thus, as the modeling period

ARTICLE IN PRESS M. Nohutcu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 72 (2010) 617–624

15

Latitude (deg)

42

10 40

38

5

36

15

42

40

623

10

38

5

36 25

30

35

40

0 45 TECU 15

42

25

35

40

0 45 TECU 15

42

40

10 40

38

5

36

30

10

38

5

36 25

30

35

40

0 45 TECU 25 Longitude (deg)

30

35

40

0 45 TECU

Fig. 5. 2-D VTEC model results for 26.09.2007, 15:30:00 (UT). VTEC maps are obtained for levels J1 ¼ 2 and J2 ¼2 for the modeling periods (UT) of: (a) 15:22:30–15:37:30 (15 min). RMSE¼ 1.0886 TECU; (b) 15:15:00–15:45:00 (30 min). RMSE ¼1.1168 TECU; (c) 15:00:00–16:00:00 (60 min). RMSE¼ 1.1433 TECU; (d) 14:30:00–16:30:00 (120 min). RMSE ¼1.1708 TECU.

5. Conclusions We presented two approaches in order to model VTEC with B-spline functions. The 2-D approach assumes the ionosphere static in a Sun-fixed reference frame for the modeling period and depends on geographic latitude and Sun-fixed longitude. The 3-D approach is presented in an Earth-fixed reference frame and

15 42 10

40 38

5

36 25

30

35

40

45 TECU

0 15

42 Latitude (deg)

gets longer, the deviations from this average state increases. Another interesting fact about the figures is that, as the modeling period gets narrower, the 2-D VTEC plots resemble the 3-D VTEC plots more and more (compare Figs. 4 and 5), which indicates that the 3-D modeling localizes the temporal variations in the ionosphere more appropriately. On the other hand, narrowing the modeling period reduces the amount of available observations for the model. B-spline solutions were compared with the solutions of the Bernese GPS Software v5.0. The Bernese GPS Software uses spherical harmonic expansion for regional or global modeling of VTEC in a Sun-fixed reference frame (Dach et al., 2007). The same observation files of the B-spline model, thus the same dataset were used within the Bernese GPS Software in order to model the VTEC depending on the Sun-fixed longitude and geographical latitude (2-D). Bernese solutions were performed for three sets of maximum degree (nmax) and maximum order (mmax) for spherical harmonics, as nmax ¼4, mmax ¼3 with 23 coefficients; nmax ¼6, mmax ¼6 with 49 coefficients; and nmax ¼12, mmax ¼8 with 149 coefficients. The corresponding RMSE values for the residuals are 1.1637 TECU, 1.1624 and 1.1625 TECU, respectively. The VTEC maps for the Bernese solutions are given in Fig. 6. Note that the Bernese GPS Software only provides the unknown coefficients of spherical harmonic expansion (not the maps) and the figures were formed with these coefficients with an external subroutine that was coded for this purpose. If we consider the RMSE values and the VTEC plots given in Figs. 3 and 6, it can be concluded that the VTEC solutions of the Bernese GPS Software are comparable with the reference solution of the 2-D B-spline solution (panel a in Fig. 3). However, on the contrary to the B-Spline solutions, increasing the maximum degree and order of spherical harmonics does not improve the detail level in VTEC maps, which is probably due to the global nature of the spherical harmonic functions.

10

40 38

5

36 25

30

35

40

45 TECU

0

15 42 10

40 38

5

36 25

30

35 40 Longitude (deg)

45 TECU

0

Fig. 6. VTEC maps which are formed with the coefficients of the Bernese GPS Software for 26.09.2007, 15:30:00 (UT). (a) nmax ¼ 4 (maximum degree), mmax ¼ 3 (maximum order), with 23 coefficients. RMSE ¼1.1637 TECU; (b) nmax ¼ 6, mmax ¼ 6 with 49 coefficients. RMSE¼ 1.1624 TECU; (c) nmax ¼ 12, mmax ¼ 8 with 149 coefficients. RMSE ¼1.1625 TECU.

considers the time dependency of the ionosphere by including a third variable for time. VTEC is split into two parts as a reference term and a correction term, where the reference term is calculated with the first level solutions of the models. For parameter estimation a robust regression algorithm, namely IRLS with a bi-square weighting function is applied in order to reduce the effects of outliers. Another iterative method, i.e. the LSQR method, is performed for the solution of linear systems. LSQR provides regularization for ill-conditioned problems and can solve very large problems with its iterative structure. B-spline approaches are applied to a real data set, which is obtained from the observations of ground-based GPS stations over Turkey. Both the RMSE values and VTEC maps indicate that the

ARTICLE IN PRESS 624

M. Nohutcu et al. / Journal of Atmospheric and Solar-Terrestrial Physics 72 (2010) 617–624

3-D solutions represent the temporal change in the ionosphere more successfully than the 2-D solutions. On the other hand, the 2-D solutions are numerically more stable as they include far less coefficients as compared to the 3-D solutions. Our B-spline solutions are compared with the solutions of the Bernese GPS Software v5.0. The VTEC solutions of the Bernese GPS Software are comparable with the reference solution of the 2-D B-spline solution. Since the presented approaches are based on Euclidean theory, their implementation is restricted to local and regional areas. However, they can be extended to larger areas or even to global scale if the B-spline functions are replaced by spherical base functions. The B-spline approach can be extended to 4-D problems (including height) to model electron density of the ionosphere in a similar manner if data are available detecting the horizontal structure of the ionosphere.

Acknowledgements This study was supported by the Scientific and Technological ¨ _ITAK) (Grant no. C Research Council of Turkey (TUB - AYDAG ¨ _ITAK-MAM for the GPS data 106Y182). The authors thank TUB provided. References Aster, R.C., Borchers, B., Thurber, C.H., 2005. Parameter Estimation and Inverse Problems. Elsevier Academic Press, USA. Bilitza, D., 2001. International Reference Ionosphere 2000. Radio Science 36 (2), 261–275. Brunini, C., Meza, A., Azpilicueta, F., Diaz, A., Van Zele, M.A., 2004. A new ionosphere monitoring technology based on GPS. Astrophysics and Space Science 290, 415–429. Chambodut, A., Panet, I., Mandea, M., Diament, M., Holschneider, M., Jamet, O., 2005. Wavelet frames: an alternative to spherical harmonic representation of potential fields. Geophysical Journal International 163, 875–899. Ciraolo, L., Azpilicueta, F., Brunini, C., Meza, A., Radicella, S.M., 2007. Calibration errors on experimental Slant Total Electron Content (TEC) determined with GPS. Journal of Geodesy 81 (2), 111–120. Cummins, D.J., Andrews, C.W., 1995. Iteratively reweighted partial least squares: a performance analysis by Monte Carlo simulation. Journal of Chemometrics 9, 489–507. Dach, R., Hugentobler, U., Fridez, P., Meindl, M. (Eds.), 2007. Bernese GPS Software Version 5.0. Astronomical Institute, University of Bern, Switzerland. DuMouchel, W.H., O’Brien, F.L., 1989. Integrating a robust option into a multiple regression computing environment. In: Proceedings of the 21st Symposium on the Interface: Computer Science and Statistics, 9–12 April 1989, Alexandria, VA, pp. 297–301. Gao, Y., Heroux, P., Kouba, J., 1994. Estimation of GPS receiver and satellite L1/L2 signal delay biases using data from CACS. In: Proceedings of KIS-94, August 30–September 2 1994, Banff, Canada, pp. 109–117.

Hansen, P.C., 1994. Regularization tools: a Matlab package for analysis and solution of discrete ill-posed problems. Numerical Algorithms 6, 1–35. Hernandez-Pajares, M., Juan, J.M., Sanz, J., 1999. New approaches in global ionospheric determination using ground GPS data. Journal of Atmospheric and Solar Terrestrial Physics 61, 1237–1247. Koch, K.R., 1999. Parameter Estimation and Hypothesis Testing in Linear Models 2nd edn. Springer, Berlin–Heidelberg–New York. Koch, K.R., Kusche, J., 2002. Regularization of geopotential determination from satellite data by variance components. Journal of Geodesy 76 (5), 259–268. Liu, Z., Gao, Y., 2003. Ionospheric TEC predictions over a local area GPS reference network. GPS Solutions 8 (1), 23–29. Lyche, T., Schumaker, L.L., 2000. A multiresolution tensor spline method for fitting functions on the sphere. SIAM Journal of Scientific Computing 22 (2), 724–746. Mannucci, A.J., Wilson, B.D., Yuan, D.N., Ho, C.H., Lindqwister, U.J., Runge, T.F., 1998. A global mapping technique for GPS-derived ionospheric total electron content measurements. Radio Science 33, 565–582. Mautz, R., Ping, J., Heki, K., Schaffrin, B., Shum, C.K., Potts, L., 2005. Efficient spatial and temporal representations of global ionosphere maps over Japan using B-spline wavelets. Journal of Geodesy 78, 660–667. Nohutcu, M., Karslioglu, M.O., Gucluer, B., Schmidt, M., Zeilhofer, C., Zhang, Z., Ergintav, S., 2008. Local modeling of VTEC using GPS observations. In: Proceedings of the TUJK Annual Scientific Meeting 2007: Monitoring and modeling of the ionosphere and troposphere, 14–16 November 2007, Ankara, Turkey, pp. 33–37. Paige, C.C., Saunders, M.A., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software 8 (1), 43–71. Schaer, S., Beutler, G., Mervart, L., Rothacher, M., Wild, U., 1995. Global and regional ionosphere models using the GPS double difference phase observable. In: Gendt G., Dick, G. (Eds.), Proceedings of the IGS Workshop 1995 on Special Topics and New Directions, 15–18 May 1995, GFZ, Potsdam, Germany, pp. 77–92. Schaer S., 1999. Mapping and predicting the Earth’s ionosphere using the Global Positioning System. Ph.D. Thesis, Astronomical Institute, University of Berne, Switzerland. Schmidt, M., 2007. Wavelet modeling in support of IRI. Advances in Space Research 39, 932–940. Schmidt, M., Fengler, M., Mayer-Gurr, T., Eicker, A., Kusche, J., Sanchez, L., Han, S.-C., 2007a. Regional gravity modeling in terms of spherical base functions. Journal of Geodesy 81, 17–38. Schmidt, M., Bilitza, D., Shum, C.K., Zeilhofer, C., 2007b. Regional 4-D modeling of the ionospheric electron content. Advances in Space Research 42, 782–790. Schmidt, M., Karslioglu, M.O., Zeilhofer, C., 2008. Regional multi-dimensional modeling of the ionosphere from satellite data. In: Proceedings of the TUJK Annual Scientific Meeting 2007: Monitoring and modeling of the ionosphere and troposphere, 14–16 November 2007, Ankara, Turkey, pp. 88–92. Stollnitz, E.J., DeRose, T.D., Salesin, D.H., 1995. Wavelets for computer graphics: a primer. Part I. IEEE Computer Graphics and Applications 15 (3), 76–84 ; Part II, IEEE Computer Graphics and Applications 15 (4), 75–85. Wielgosz, P., Grejner-Brzezinska, D., Kashani, I., 2003. Regional ionosphere mapping with kriging and multiquadratic methods. Journal of Global Positioning System 2 (1), 48–55. Yao, Z.S., Roberts, R.G., Tryggvason, A., 1999. Calculating resolution and covariance matrices for seismic tomography with the LSQR method. Geophysical Journal International 138, 886–894. Zeilhofer, C., 2008. Multi-Dimensional B-Spline Modeling of Spatio-Temporal Ionospheric Signals. German Geodetic Commission, Series A, 123, Munich. Zeilhofer, C., Schmidt, M., Bilitza, D., Shum, C.K., 2009. Regional 4-D modeling of the ionospheric electron density from satellite data and IRI. Advances in Space Research 43, 1669–1675.