84 Chapter 5. Modeling the ionosphere considered to simulate a realistic geometrical ray path distribution. For one hour data between 06:00 UT and 07:00 UT at 1 July 2002, the geometry of five GPS satellite passes is shown by Fig. 5.3.
Satellite passes are depicted as blue lines and receivers as red triangles where pass start- and end-points are connected by dashed lines. From these stations, totally 2,404 synthetic TEC observations according to the epochs provided by their RINEX data are generated with 30 s sampling according to Eq. (5.3). The Chapman parameters are derived from IRI data while the integration is adjusted with N = 9 and 1 km step width fixed for all three layers to obtain a synthetic reference solution for the validation of the internal accuracy. The Root Mean Square (RMS) of the reference solution for this selected testbed yields 21.85 TECU.
Afterwards, integrations are performed with varying vertical step widths of 20 km, 40 km, 60 km and 80 km and different quadrature order N ∈ {3,5,7,9}. Figure 5.4 depicts the comparison of the RMS values related to the TEC differences (green bars) and the processing time (blue bars). The darker the color, the smaller the step width and consequently, the computation time decreases with decreasing quadrature order and increased step width. The RMS on the contrary increases with the step width and decreasing order. In this internal validation, it can be seen that the RMS values are permanently small while the processing time must be considered as the limiting factor, in particular since the results here are based on only 2,404 observations. Typically, a maximum quadrature level of N = 6 was considered in this work. A vertical step width of 20 km was assigned for the second layer to account for increased ionospheric gradient and larger step widths between 60 km and 80 km are commonly chosen for the first and third layer.
Figure 5.4: Comparison of processing effort in computation time (blue bars) with the RMS of TEC differences (green bars) for different quadrature orders and integration step widths in a simulated scenario.
5.3. Parameter estimation 85 1999; Niemeier, 2008). The linear deterministic model describing the functional relation between the observations and unknown target parameters is defined by
(N×1)l + e
(N×1)= A
(N×U) x
(U×1) (5.7)
where l = [l1,l2, . . . ,lN]T contains the input observations, e = [e1,e2, . . . ,eN]T denotes randomly distributed measurement errors andAis the real-numbered coefficient or design matrix. The unknown target parameters are collected in the vectorx= [x1,x2, . . . ,xU]T. NandUare introduced here as the number of observations and unknowns, respectively. In case that the dependencies f(x) = l +e are given as non-linear functions, a Taylor39series expansion must be applied to linearize each equation
fi around the pointx0of initial values with fi(x0+∆x) = fi(x0)+
∂fi(x)
∂x
T
∆x (5.8)
under negligence of second O (∆x)2
and higher order terms. The partial derivatives are collected in the Jacobian40matrix A.
Associated with Eq. (5.7), the stochastic part of the Gauß-Markov model is defined as Σll
(N×N) = σ02 Pll−1
(N×N)
= σ20 Qll
(N×N)
(5.9) where Pll is the positive-definite observation weight matrix and σ02 the initially unknown variance factor. Σll identifies the covariance matrix that contains the standard deviations σi and correlation coefficients ρof the observations such that
Σll =
σ12 ρ12σ1σ2 · · · ρ1Nσ1σN
ρ21σ2σ1 σ22 ...
... . . .
ρN1σNσ1 · · · σ2N
. (5.10)
The inverse ofPll is often denoted as cofactor matrixQll. In case that no correlations between obser-vations are taken into account, the matrices Pll,Qll and consequently alsoΣll have a diagonal struc-ture including only the individual observation variances along the diagonal. For equally weighted, uncorrelated observations it follows thatPll can be set to the unit matrixI.
In Section 5.2, the fundamental observation equations for the electron density and TEC have been introduced by means of the Eqs. (5.2) and (5.3). Due to the exponential term contained in the Chapman function, the dependencies are non-linear and need to be linearized to comply with the conditions of the Gauß-Markov model. According to the system of linear relations of Eq. (5.7), each measurement can be described by means of Eq. (5.8) with
li+ei = fi(x0+∆x) = fi(x0)+
U
X
j=1
∂fi
∂xj
!
x0
∆xj (5.11)
wherex0indicates the linearization point that is defined by a priori information. Now, a linear system of equations with
L+e= A∆x (5.12)
39Brook Taylor (18.08.1685-29.12.1731) was a British mathematician and member of the royal society.
40After Carl Gustav Jacob Jacobi (10.12.1804-18.02.1851), a German mathematician.
86 Chapter 5. Modeling the ionosphere can be established including the reduced observation vectorL =l− f(x0) and the design matrix
(N×U)A =
a11 a12 · · · a1U a21 . . . ...
...
aN1 · · · aNU
where ai j = ∂fi
∂xj
!
x0
. (5.13)
Following the principle of minimum as postulated by the method of least squares, the errors
e= A∆x−L shall satisfy the minimum condition for eTPlle. (5.14) A minimum can be found at the zeros of the first derivative ofeTPlle which is true for ATPllA∆Dx =! ATPllL. The product ATPllAis commonly identified as normal equation matrix, henceforth denoted as NEQ. IfNEQ is of full rank and thus invertible, a solution for the unknown parameters can be found with
∆Dx =
ATPllA−1
ATPllL (5.15)
as the Best Linear Unbiased Estimation (BLUE) (Koch, 1999).
Observation equation
The theory shall now be applied to determine the electron density or TEC, respectively. The emphasis is on the modeling of the F2 layer by estimating the F2 Chapman parameters. The Eqs. (5.2) and (5.3) describe the relation between the Chapman parameters and the observations. The plasmasphere parameters have been approximated by constant values with physical meaning. Following Jakowski (2005), a reasonable value of 104km for HP (in the case of h > hmF2, otherwise 10 km) and the assumption of proportionality betweenN0P andN mF2 can be taken into account. The following steps can of course also be applied for other ionospheric layers and the parameter basis may be extended by additional unknowns.
With
fi(x) = Ne(N mF2,hmF2,HF2) or fi(x) =
xs
Z
xr
Ne(N mF2,hmF2,HF2)ds, (5.16) the observation vectorl contains electron density or TEC measurements depending on the input data of the ionospheric sounding technique.
Under consideration of the linearization pointx0 = (N mF2|0,hmF2|0,HF2|0)T, the electron density can be modeled as
li+ei−Ne,i(x0)
| {z }
Li
=
∂Ne,i
∂N mF2 x0
∆N mF2+
∂Ne,i
∂hmF2 x0
∆hmF2+
∂Ne,i
∂HF2 x0
∆HF2 (5.17) taking the Chapman profile function into account to represent the vertical distribution. For modeling the horizontal and temporal variation, a B-spline parametrization as specified in Section 2.3 is applied.
Polynomial B-splines have, for instance, been considered by Schmidt et al. (2008) and Zeilhofer et al.
(2008) for modeling the regional GPS TEC. P-I: Limberger et al. (2013) incorporated the Chapman profile into the approach to represent regional electron density from occultation data. Furthermore, Schmidt et al. (2011) and P-IV: Dettmering et al. (2014b) applied a combination of trigonometric and polynomial B-splines for modeling the global TEC distribution based on GPS and DORIS measure-ments, respectively.
5.3. Parameter estimation 87 Each Chapman key parameter is decomposed into tensor products of three one-dimensional B-splines with respect to latitude ϕ, longitude λ and time t. The reference systems are typically defined in an Earth- or Sun-fixed frame and therefore ϕidentifies the geographical or geomagnetic latitude and λ the geographical or Sun-fixed longitude (= local time). Based on the general notation κh with h∈ {1,2, . . . ,H}which is henceforth used to identify the key parameters, the B-spline representation yields
κh(λ, ϕ,t) =
K1
P
k1=1 K2
P
k2=1 K3
P
k3=1dkJ1,J2,J3
1,k2,k3 φkJ1
1(ϕ) φJk2
2(λ)φJk3
3(t) regional
K1
P
k1=1 K2
P
k2=1 K3
P
k3=1
dkJ1,J2,J3
1,k2,k3 φkJ1
1(ϕ)TkJ2
2(λ) φkJ3
3(t) global
. (5.18)
For regional applications, the F2 Chapman parameters, N mF2, hmF2 andHF2, are thus represented as
κ1(ϕ, λ,t)R = N mF2(ϕ, λ,t)R =
K1
X
k1=1 K2
X
k2=1 K3
X
k3=1
dJk1,J2,J3
1,k2,k3
N mF2 φkJ1
1(λ) φJk2
2(ϕ)φJk3
3(t) (5.19) κ2(ϕ, λ,t)R = hmF2(ϕ, λ,t)R =
K1
X
k1=1 K2
X
k2=1 K3
X
k3=1
dJk1,J2,J3
1,k2,k3
hmF2 φJk1
1(λ) φkJ2
2(ϕ) φkJ3
3(t) (5.20) κ3(ϕ, λ,t)R = HF2(ϕ, λ,t)R =
K1
X
k1=1 K2
X
k2=1 K3
X
k3=1
dJk1,J2,J3
1,k2,k3
HF2 φJk1
1(λ)φJk2
2(ϕ) φkJ3
3(t) (5.21) wheredkJ1,J2,J3
1,k2,k3 are the unknown scaling coefficients that need to be determined. For global modeling, φkJ2
2(ϕ) may be replaced by trigonometric B-splinesTkJ2
2(ϕ). Further details about the application of B-splines for ionosphere modeling can be found in Section 2.3. It shall however be noted once more, that the number of unknowns depends on the numberKof spline functions and therefore on the model resolution defined by the level J. It follows that the unknowns collected inx areU = QJ
j=1Kj B-spline series coefficients.
Under consideration of the same B-spline levels for the longitude and time (J2= J3= 2) but a diff er-ent level for the latitude (J1 = 3), the number of unknowns for the regional model with polynomial B-splines yields, for instance,U =360 coefficients per key parameter.
To account for B-spline coefficients as target parameters, Eq. (5.17) needs to be adapted. The modi-fication leads to
li+ei−Ne,i(κ(x0))
| {z }
Li
=
H
X
h=1
* ,
∂Ne
∂κh
κ(x0)
∂κh
∂dκh x0
+
-T
i
∆dκh (5.22)
where x0 = [d0,κ1,d0,κ2, . . . ,d0,κH]T contains initial values of the series coefficients related to the Chapman key parameters. All key parameters are contained in the vectorκ.
The model equation for the integrated electron density along the ray path between transmitter and receiver follows instantly from Eq. (5.3) and according to Eq. (5.22) can be obtained as
li+ei−
xs
Z
xr
Ne,i(κ(x0))ds
| {z }
Li
=
xs
Z
xr
* . ,
H
X
h=1
* ,
∂Ne
∂κh
κ(x0)
∂κh
∂dκh
x0
+
-T
i
∆dκh+ /
-ds. (5.23)
88 Chapter 5. Modeling the ionosphere
Partial derivatives of F2 model parameters
The partial derivatives of Eq. (5.2) with respect to the F2 Chapman parameters are part of the right-hand sides in Eq. (5.22) and Eq. (5.23) and will be placed in the design matrix A. Therefore, the corresponding terms for N mF2,hmF2 and HF2 are summarized in the following.
1. Partial derivatives of theα-Chapman function for N mF2 and the B-spline coefficients:
∂Ne,i
∂N mF2 =exp 1
2 1− h−hmF2
HF2 −exp −h−hmF2 HF2
! ! !
(5.24)
∂Ne,i
∂N mF2
∂N mF2
∂
dkJ1,J2,J3
1,k2,k3
N mF2
= ∂Ne,i
∂N mF2 ·φJk1
1(λ) φkJ3
3(t)
φkJ2
2(ϕ) regional TkJ2
2(ϕ) global (5.25)
2. Partial derivatives of theα-Chapman function forhmF2 and the B-spline coefficients:
∂Ne,i
∂hmF2 =
N mF2 2HF2
∂Ne,i
∂N mF2
1−exp
−h−hmF2HF2 + N0PHP ∂N0P∂Ne,i ifh ≥ hmF2
N mF2 2HF2
∂Ne,i
∂N mF2
1−exp
−h−hmF2HF2 − N0PHP ∂N0P∂Ne,i else (5.26)
∂Ne,i
∂hmF2
∂hmF2
∂
dkJ1,J2,J3
1,k2,k3
hmF2
= ∂Ne,i
∂hmF2φkJ1
1(λ) φJk3
3(t)
φJk2
2(ϕ) regional TkJ2
2(ϕ) global (5.27)
3. Partial derivatives of theα-Chapman function for HF2 and the B-spline coefficients:
∂Ne,i
∂HF2 = N mF2 ∂Ne,i
∂N mF2
h−hmF2
2(HF2)2 1−exp −h−hmF2 HF2
! !
(5.28)
∂Ne,i
∂HF2
∂HF2
∂
dkJ1,J2,J3
1,k2,k3
HF2
= ∂Ne,i
∂HF2φkJ1
1(λ) φkJ3
3(t)
φkJ2
2(ϕ) regional TkJ2
2(ϕ) global (5.29)
The functional dependencies between measurements and model parameters define the basis of the ad-justment system including the unknown parameters that will be estimated. Therefore, the observation equations for those techniques presented in the Chapters 3 and 4 are summarized in the following and the structure of the associated design matrix is provided. It should be noticed, that the unknowns introduced here are only the very essential ones to be considered for regional modeling. For global modeling, additional B-spline constraints as described in Section 2.4 may be taken into account. Be-sides, the estimation of bias parameters with respect to a reference technique are recommended to absorb relative biases between different missions and observation techniques.
Ionospheric GNSS radio occultations
The electron density retrieval procedure described in Chapter 4 allows to obtain electron density profiles including discrete Ne values derived from STEC observations between a GNSS transmitter
5.3. Parameter estimation 89 and a LEO receiver. After a relatively complex preprocessing and data retrieval, the retrieved electron density values can be modeled by
Ne+e= N mF2 exp 0.5 1− h−hmF2
HF2 −exp −h−hmF2 HF2
! ! !
+N0P −|h−hmF2|
HP
! (5.30) resulting directly from Eq. (5.2) under consideration of the extendedα−Chapman layer. As aforemen-tioned, the design matrixAcontains the partial derivatives and the vector of unknowns the corrections of the B-spline series coefficients. It follows that
AIRO=
∂Ne,1
∂κ1 ∂κ1
∂d · · · ∂N∂κe,1
H
∂κH
∂Ne,2 ∂d
∂κ1 ∂κ1
∂d · · · ∂N∂κe,2
H
∂κH
... . . . ... ∂d
∂Ne,N
∂κ1 ∂κ1
∂d · · · ∂N∂κe,N
H
∂κH
∂d
(5.31)
and
∆DxIRO = f
∆dDκ1,1,∆dDκ1,2, . . . ,∆dDκ2,1,∆dDκ2,2, . . . ,∆dDκH,1,∆dDκH,2, . . .gT
. (5.32)
GNSS
In case of GNSS, the integrated electron density along a slant ray path as described in Section 3.2 is observed by means of leveled carrier phase observations, i.e., Eq. (5.3) must be taken into account where additional instrumental biases, the DCBs, for the satellites and receivers need to be considered.
The observation equation yields STEC+e =
xs
Z
xr
Neds+DCBr +DCBs
=
xs
Z
xr
N mF2 exp 0.5 1− h−hmF2
HF2 −exp −h−hmF2 HF2
! ! !
+N0P −|h−hmF2|
HP
! ds
+DCBr +DCBs. (5.33)
It is assumed, that an approximation for the carrier phase bias has already been determined by means of a suitable leveling procedure. With the DCBs as additional unknowns, the system blows up where Aand∆Dxare expanded to
AGNSS=
∂Ne,1
∂κ1 ∂κ1
∂d · · · ∂N∂κe,1
H
∂κH
∂d ∂Ne
∂DCBr1 · · · ∂DCB∂Ne
r Nr
∂Ne
∂DCBs1 · · · ∂DCB∂Nes Ns
∂Ne,2
∂κ1 ∂κ1
∂d · · · ∂N∂κe,2
H
∂κH
∂d ∂Ne
∂DCBr1 · · · ∂DCB∂Ne
r Nr
∂Ne
∂DCBs1 · · · ∂DCB∂Nes Ns
... . . . ... ... . . . ... ... . . . ...
∂Ne,N
∂κ1 ∂κ1
∂d · · · ∂N∂κe,N
H
∂κH
∂d ∂Ne
∂DCBr1 · · · ∂DCB∂Ne
r Nr
∂Ne
∂DCBs1 · · · ∂DCB∂Nes Ns
(5.34) and
∆DxGNSS =
"
∆dD dcbM
#
(5.35)
90 Chapter 5. Modeling the ionosphere with
∆dD= f
∆dDκ1,1,∆dDκ1,2, . . . ,∆dDκ2,1,∆dDκ2,2, . . . ,∆dDκH,1,∆dDκH,2, . . .gT
(5.36) dcbM =
DCB[r1,· · · ,DCB[rNr,DCB[s1, . . . ,DCB[sNs T
. (5.37)
NrandNsindicate the number of distinct receivers and satellites considered in the system. The partial derivatives for the DCBs of a specific observation are one for the satellite-receiver pair and otherwise zero. For the reason, that all DCBs may be shifted by a common offset, a zero-mean assumption for relative constraining is commonly introduced for the satellite DCBs with
Ns
X
i=s1
DCBi
=! 0. (5.38)
This constraint can, for instance, be considered by means of an additional observation equation in-cluded at the normal equations level.
DORIS
The derivation of ionospheric signals from DORIS measurements has been described in Section 3.3.
A strong drawback of DORIS data is the high noise level in the pseudorange observations causing that standard procedures for the carrier phase leveling cannot be applied. Thus, the observation equation
STEC+e =
xs
Z
xr
Neds+CPB
=
xs
Z
xr
N mF2 exp 0.5 1− h−hmF2
HF2 −exp −h−hmF2 HF2
! ! !
+N0P −|h−hmF2|
HP
!
ds+CPB (5.39) is related to unleveled data and still contains a CPB term including, e.g., the PWU and in particular the phase ambiguities. As a preliminary strategy to account for the CPB, the TEC data of IGS GIMs have been exploited by Dettmering et al. (2014a); P-IV: Dettmering et al. (2014b) as a reference to level the DORIS STEC. The entries of Aandxˆ are basically related to the series coefficients as shown in Eq. (5.31) and Eq. (5.32). However, the consideration of an additional offset with respect to other techniques and missions is recommended to account for systematic biases.
Altimetry
From altimetry, vertical range errors are provided that can directly be converted into VTEC. Thus, the observation equation reads
VTEC+e =
xs
Z
xr
Nedh
=
xs
Z
xr
N mF2 exp 0.5 1− h−hmF2
HF2 −exp −h−hmF2 HF2
! ! !
+N0P −|h−hmF2|
HP
! dh.
(5.40) As for occultations and DORIS, the minimum unknowns are defined with the B-spline coefficients.
The existence of biases with few TECU between altimetry and GNSS is well known. The relative leveling between the techniques and also altimetry missions is therefore highly recommended.