Computing unambiguous TEC and ionospheric delays using only carrier
phase data from NOAA's CORS network

Dru A. Smith

National Geodetic Survey, NOS/NOAA

Silver Spring, Maryland 20910, USA.

***Abstract*-A new method for computing absolute (unambiguous) levels of
Total Electron Content (TEC) and subsequently the L1 and L2 phase delays
of GPS is presented. Unlike previous computational methods, this one
relies solely upon dual frequency, ambiguous carrier phase data without
any reliance on pseudo-range, a-priori values or other external
information. The only requirements for this method are that the
ionosphere is assumed to lie in a two-dimensional shell of constant
ellipsoidal height, and that the GPS data come from a network of ground
stations, geographically separated so as to allow satellites to be
viewed by a variety of stations at overlapping times.**

**The usefulness of this method and its application toward nowcasting
and forecasting of the ionosphere are also discussed.**

[I. Introduction]{.smallcaps}

It is a well accepted fact that dual-frequency carrier phase data is
currently the most accurate data for computing precise positions with
the Global Positioning System. However, both integer ambiguities and the
Total Electron Content (TEC) along the receiver-satellite path are
frequency dependent; and both of these values prevent an absolute and
immediate computation of precise position from carrier phase data alone.

Many methods exist for either removing the ionosphere, or mixing in
pseudo-range data, or computing both ambiguities and ionosphere delay
using lengthy observing sessions. However, in this paper, a method for
unambiguously computing the absolute level of TEC, from a network of
ground stations and only using ambiguous carrier phase data is being
presented. With TEC data computed, the burden of determining integer
ambiguities should be significantly reduced \[1,2\]. In most studies,
however, the presumption is that data, including the ionosphere, are
processed in double difference mode. This study goes beyond that,
actually computing absolute values of TEC from ambiguous carrier phase
data.

While the method presented does require the assumption of a
two-dimensional "shell model" of the ionosphere, independent checks
confirm that the data coming from this method are generally accurate,
though with some study still needed to resolve some discrepancies with
coarse global ionosphere models. The method presented is fast, and can
be kept updated over time using a sequential least squares adjustment as
the data change (as new satellites rise and old ones set).

[II. Assumptions and Definitions]{.smallcaps}

This paper makes use of the simple shell model of the ionosphere. In
this model, the three-dimensional ionosphere, with TEC measured in
electrons/m^3^, is compressed to a shell at a fixed ellipsoidal height
(often picked as the location approximating the largest real density of
electrons). The compression takes place in the direction perpendicular
to the shell and leaves the ionosphere with a TEC that is expressed in
electrons/m^2^ (see Figure 1). For this paper, the shell is fixed at a
height of 300 km.

![](media/image1.png){width="2.875in" height="2.196527777777778in"}

Fig. 1. Sketch of the relationship between 3-D electron density and the
compression to a 2-D shell, showing the definition of a pierce point.

With the ionosphere shell model in place, consider next the line
connecting a GPS receiver with a GPS satellite, at any given data
collection epoch. The intersection of this line with the ionosphere will
be called a [pierce point]{.underline} (see Figure 1). If one maps the
pierce points generated this way, from the rise of satellite above the
horizon to its setting below the horizon, a [track]{.underline} is
generated, such as that seen in Figure 2.

![](media/image2.png){width="3.370833333333333in"
height="3.160416666666667in"}

Fig. 2. Track 10610 formed when station SBT1 tracks GPS SV 3 from rise
to set. Day 193 (July 12) of 2002. Station SBT1 is shown by the large
asterisk.

Consider now a close-up on the ionosphere, at a particular pierce point
(see Figure 3).

![](media/image3.png){width="3.1458333333333335in" height="2.53125in"}

Fig. 3. Close-up of the ionosphere shell at a pierce point, showing the
difference between the two electron densities TECS and TECR

$$An important distinction will be made between the density of electrons
seen along the receiver-satellite direction (TECS) and the density of
electrons seen in a direction radial to the ionosphere shell (TECR).
They are related by the piercing angle, z' as:

$\text{TECS} = \text{TECR}\text{/cos}z'$ (1)

It has long been understood that the phase advance (or negative delay)
experienced by GPS receivers is a function of the number of electrons
through which the L1 and L2 carrier waves pass. That is, the number of
cycles of advance, due to the ionosphere, can be written as a function
of TECS and frequency as:

$I_{1} = - \frac{\kappa}{f_{1}c}\text{TECS} = (0\text{.}\text{853})\text{TECS}$
(2)

> $I_{2} = - \frac{\kappa}{f_{2}c}\text{TECS} = (1\text{.}\text{095})\text{TECS}$
> (3)

where:

> I~1~,I~2~ = cycles of phase delay in L1 and L2 respectively

c = speed of light \[299,792,458 m / s\]

> f~1~, f~2~ = frequency of L1 and L2 respectively \[1,557,420,000 and
> 1,227,600,000 cyc / s\]
>
> TECS = density of electrons in receiver-satellite direction in TECU
> (10^16^ elec / m^2^)

κ = 40.3 x 10^16^ (m cyc^2^ / TECU / s^2^)

While the ionospheric delay experienced at a receiver is a function of
TECS, it would be most helpful to somehow generate a map of TECR, so
that any vector connecting a GPS satellite to a GPS receiver, piercing
the ionosphere, can have its TECS value (and thus its ionospheric delay)
computed by simply knowing TECR and computing the piercing angle z' from
the satellite orbit.

As a final definition, consider the case when two tracks come "near" one
another. In the field of altimetry, if they actually crossed
geographically this would be called a "crossover". In this study,
however, the term [crossover]{.underline} will not be applied to strict
crossing, but rather to any two pierce points, on separate tracks, that
fall within some acceptably small tolerance (in *both* space *and* time)
of one another.

[III. The Mathematical Model]{.smallcaps}

At any epoch "i", the relationship between the geometric range (^i^d)
and the measured range (^i^r) from a GPS receiver to satellite can be
expressed (for L1 and L2 respectively) as:

$r_{1} = d + c\Delta t + \frac{I_{1}}{\lambda_{1}} + T + m_{1}$ (4)

$r_{2} = d + c\Delta t + \frac{I_{2}}{\lambda_{2}} + T + m_{2}$ (5)

where

> r~1~, r~2~ = Measured range to satellite on L1 and L2 respectively

d = Geometric range to satellite

> c∆t = Range error due to clock errors
>
> λ~1~, λ~2~ = Wavelengths of L1 and L2 respectively
>
> T = Range error due to troposphere
>
> m~1~,m~2~ = Range error due to carrier phase multipath error on L2 and
> L2 respectively

Now, presume that a GPS receiver is tracking carrier phase data on L1
and L2. At acquisition of L1 and L2, it will compare its internal
oscillator with the number of L1 and L2 cycles received and report a
more or less nonsensical (i.e. "ambiguous") number of cycles. Presuming
it does not lose lock, it will continue to generate an internal L1and L2
cycle count, and compare the changes in cycles of its internal count to
those cycles received from the GPS satellite. The difference between
received and internal cycles can be interpreted as the change in
measured range (r) due to satellite motion. Equivalently it could be
considered the change in geometric range (d) plus the change to all
error sources. Thus between two epochs i and j:

$\begin{array}{r}
{\Delta\phi}_{1}\lambda_{1} = {\Delta r}_{1} = \\
\Delta d + c\text{ΔΔ}t + \frac{{\Delta I}_{1}}{\lambda_{1}} + \Delta T + {\Delta m}_{1}
\end{array}$ (6)

$\begin{array}{r}
{\Delta\phi}_{2}\lambda_{2} = {\Delta r}_{2} = \\
\Delta d + c\text{ΔΔ}t + \frac{{\Delta I}_{2}}{\lambda_{2}} + \Delta T + {\Delta m}_{2}
\end{array}$ (7)

where ^i,j^∆φ~1~ and ^i,j^∆φ~2~are the differences, in cycles, between
the number of L1 or L2 cycles generated inside the receiver with the
number of cycles received from the GPS satellite, over the time period
from epoch i to epoch j. (This is basically the difference between any
two observed values of L1 or L2 between two epochs). Now, if one
subtracts equation 7 from equation 6 then only frequency dependent terms
will remain:

$\begin{array}{r}
{\Delta\phi}_{1}\lambda_{1} - {\Delta\phi}_{2}\lambda_{2} = \\
(\frac{{\Delta I}_{1}}{\lambda_{1}} - \frac{{\Delta I}_{2}}{\lambda_{2}}) + \text{ΔΔ}m_{1,2}
\end{array}$ (8)

Going forward, the assumption will be made that m~1~ and m~2~ are small,
and also that the difference between m~1~ and m~2~ is small enough to be
neglected \[3,4\]. Then, applying equations 1, 2 and 3 to equation 8 and
solving for the TECS term we end up with the following:

$\begin{array}{r}
\Delta\text{TECS} = \text{TECS} - \text{TECS} = \\
\left( \frac{1}{\kappa} \right)\left( \frac{1}{f_{1}^{2}} - \frac{1}{f_{2}^{2}} \right)^{- 1}\left( {\Delta\phi}_{1} - {\Delta\phi}_{2} \right)
\end{array}$ (9)

This equation shows that if we take the difference (in time) of the
difference (in frequency) of L1 and L2 cycle changes (between internal
and received) we have an unambiguous measure of the difference (in time)
of the electron density as seen along the receiver-satellite path
(TECS). What this means is that if we plot a curve of TECS versus time
for a particular track, we know the *shape* of TECS for that track, but
we are *missing a single bias* which would define the absolute values of
TECS for that track. This is exemplified in Figure 4 where five (of an
infinite number) of TECS curves have been drawn for one particular track
(designated number 10,610 of some 12,585 tracks during the day). Note
that all 5 curves have the same shape, but with an unknown absolute
value.

![](media/image4.png){width="3.3541666666666665in"
height="2.5416666666666665in"}

Fig. 4. Five possible versions of TECS for track 10610 showing identical
shape, but variable absolute values.

Remember that the goal of this paper is to model TECR values, not TECS.
An equation for ∆TECR as a function of ∆TECS can be written as:

$\begin{array}{r}
\Delta\text{TECR} = \text{TECR} - \text{TECR} = \\
\text{TECS}\text{cos}z' - \text{TECS}\text{cos}z' = \\
\Delta\text{TECS}\text{cos}z' + \text{TECS}(\text{cos}z' - \text{cos}z')
\end{array}$ (10)

The disconcerting conclusion of equation 10 is that we may not simply
determine the unambiguous temporal changes of TECR along a track the way
we did for TECS (the dependence of ^i,j^∆TECR upon ^i^TECS causes the
ambiguity). This ambiguity, arising from the cos z' mapping between TECS
and TECR prevents us from determining even the shape of TECR curves, let
alone their absolute values. This is shown in Figure 5 where the five
TECR curves corresponding to the five TECS curves of Figure 4 are
plotted. Notice that they are not parallel curves (like TECS was), but
rather have shapes that depend upon their absolute values.

While this at first may seem a point of difficulty, as it turns out the
non-uniqueness of the shape of TECR curves is exactly the key
information needed to determine unambiguous TECR data from a network of
GPS stations. The basic idea is explained below.

![](media/image5.png){width="3.3541666666666665in"
height="2.5416666666666665in"}

Fig. 5. Five possible versions of TECR for track 10610 (corresponding to
the TECS curves of figure 4) showing differing shapes dependent upon
absolute values.

Consider the case of 3 GPS tracks (call them 1,2 and 3), generated by 3
GPS/receiver combinations, and distributed in such a way that the 3
tracks have 3 crossovers (call them A, B and C), effectively forming a
"triangle" on the ionosphere shell (see Figure 6).

![](media/image6.png){width="2.875in" height="2.3854166666666665in"}

Fig. 6. A "triangle" formed by 3 tracks and 3 crossovers.

As previously stated, we know the shape of the TECS curve along each of
the three tracks, but do not know the absolute value. Put another way,
each of the 3 tracks (1, 2 and 3) has 1 unknown bias (b~1~, b~2~ and
b~3~). In fact, let's define b~i~ as the actual value of TECS on track i
at the very first point on that track. If we knew that bias, we could
compute the absolute values of TECS (because we know ∆TECS along the
track) and knowing the absolute levels of TECS we can compute the
absolute levels of TECR. So in this example, 3 tracks mean 3 unknowns.
But how do we solve for them? As it turns out, each crossover acts as a
unique constraint for the system. The three crossovers at A, B and C
provide the following constraints on the system:

$_{1}^{A}\text{TECR} =_{3}^{A}\text{TECR}( =_{}^{A}\text{TECR})$ (11)

$_{1}^{B}\text{TECR} =_{2}^{B}\text{TECR}( =_{}^{B}\text{TECR})$ (12)

$_{2}^{C}\text{TECR} =_{3}^{C}\text{TECR}( =_{}^{C}\text{TECR})$ (13)

That is, the TECR values at a crossover must be unique, independent of
track. Notice that while all of our actual track information is in the
form of ∆TECS values, the actual constraints on the system are on TECR.
So, in the notation used above, let's write out the conversion from TECR
to TECS and see how the unknown biases fit into the picture:

$\begin{array}{r}
_{1}^{A}\text{TECR} =_{1}^{A}\text{TECS}\text{cos}_{1}^{A}z' = \\
(b_{1} +_{1}^{A}\Delta\text{TECS})\text{cos}_{1}^{A}z'
\end{array}$ (14)

$\begin{array}{r}
_{3}^{A}\text{TECR} =_{3}^{A}\text{TECS}\text{cos}_{3}^{A}z' = \\
(b_{3} +_{3}^{A}\Delta\text{TECS})\text{cos}_{3}^{A}z'
\end{array}$ (15)

$\begin{array}{r}
_{1}^{B}\text{TECR} =_{1}^{B}\text{TECS}\text{cos}_{1}^{B}z' = \\
(b_{1} +_{1}^{B}\Delta\text{TECS})\text{cos}_{1}^{B}z'
\end{array}$ (16)

$\begin{array}{r}
_{2}^{B}\text{TECR} =_{2}^{B}\text{TECS}\text{cos}_{2}^{B}z' = \\
(b_{2} +_{2}^{B}\Delta\text{TECS})\text{cos}_{2}^{B}z'
\end{array}$ (17)

$\begin{array}{r}
_{2}^{C}\text{TECR} =_{2}^{C}\text{TECS}\text{cos}_{2}^{C}z' = \\
(b_{2} +_{2}^{C}\Delta\text{TECS})\text{cos}_{2}^{C}z'
\end{array}$ (18)

$\begin{array}{r}
_{3}^{C}\text{TECR} =_{3}^{C}\text{TECS}\text{cos}_{3}^{C}z' = \\
(b_{3} +_{3}^{C}\Delta\text{TECS})\text{cos}_{3}^{C}z'
\end{array}$ (19)

Substituting equations 14-19 into equations 11-13, and converting to
matrix (Y=A X) form yields:

$\begin{bmatrix}
_{1}^{A}\Delta\text{TECS}\text{cos}_{1}^{A}z' -_{3}^{A}\Delta\text{TECS}\text{cos}_{3}^{A}z' \\
_{1}^{B}\Delta\text{TECS}\text{cos}_{1}^{B}z' -_{2}^{B}\Delta\text{TECS}\text{cos}_{2}^{B}z' \\
_{2}^{C}\Delta\text{TECS}\text{cos}_{2}^{C}z' -_{3}^{C}\Delta\text{TECS}\text{cos}_{3}^{C}z'
\end{bmatrix} = \begin{bmatrix}
 - \text{cos}_{1}^{A}z' & 0 & + \text{cos}_{3}^{A}z' \\
 - \text{cos}_{1}^{B}z' & + \text{cos}_{2}^{B}z' & 0 \\
0 & - \text{cos}_{2}^{C}z' & + \text{cos}_{3}^{C}z'
\end{bmatrix}\begin{bmatrix}
b_{1} \\
b_{2} \\
b_{3}
\end{bmatrix}$ (20)

This system of linear equations has a similarity to those used, for
example, in solving altimetric crossovers or geodetic leveling problems.
In the case of altimetry and geodetic leveling however, the A matrix
contains only +1, -1 and 0 and is thus non-invertible without additional
constraints. In contrast, the cosines in our A matrix allow for an
inversion. In simplest terms, this means that even though we only know
∆TECS (that is, the shape or temporal changes of TECS) on a track and we
only know that TECR must be identical at the crossovers (independent of
track), we can compute one unique absolute value of TECS at the
beginning of each track (b~1~, b~2~ and b~3~). Knowing this value, we
are then able to compute all TECS values for the entire track and
subsequently all absolute TECR values for the entire track.

Therefore, within the confines of our opening assumptions, an
unambiguous solution for TECR, based solely on ambiguous carrier phase
data, has been found. It is entirely due to the (cos z') mapping between
TECS and TECR that this is possible. As it turns out, this mathematical
relationship holds for any number of tracks that form a closed polygon,
not just for triangles, since the number of sides (tracks) of a polygon
must equal the number of vertices (crossovers), and thus would provide a
unique and solvable system.

However, while the above triangular example offers a unique solution, it
offers no redundancy. But if the situation is made just slightly more
complex, say by adding a 4th track, crossing tracks 1 and 2 at points D
and E respectively (forming another triangle, with common vertex B) we
begin to build up redundancy. See Figure 7.

![](media/image7.png){width="2.875in" height="2.4270833333333335in"}

Fig. 7. Adjoining "triangles" showing 4 tracks (unknowns) and 5
crossovers (observations) and thus redundancy

In this case, we now have 4 tracks (with unknowns b~1~, b~2~, b~3~ and
b~4~) but 5 crossovers (observation equations), and thus a Least-Squares
Adjustment (LSA) can be performed on this redundant system. This buildup
of redundancy grows as the number of tracks grows, since each new track
can be expected to cross more and more existing tracks, adding multiple
crossovers (observations) to the system while only adding 1 unknown bias
for the new track.

Given the fact that a complex network of tracks should yield more
crossovers than tracks, it seems that the ideal application of this
model is toward a large network of dual frequency carrier phase GPS
receivers, which are continuously operating and spaced geographically so
that their tracks cross, but offer some wide spatial distribution in the
ionosphere.

[IV. The CORS Network]{.smallcaps}

The National Geodetic Survey (NGS), part of the National Oceanic and
Atmospheric Administration (NOAA), coordinates a network of continuously
operating reference stations (CORS) that provide GPS carrier phase and
code range measurements in support of 3-dimensional positioning
activities throughout the United States and its territories (see Figure
8).

![](media/image8.png){width="3.2291666666666665in"
height="2.5833333333333335in"}

Fig. 8. The CORS Network

The CORS system benefits from a multi-purpose cooperative endeavor
involving many government, academic, commercial and private
organizations. New sites are evaluated for inclusion according to
established criteria. The data are collected at various rates (1, 5, 10
or 30 seconds) and then transmitted (with varying latencies) to NGS for
coordinate computation, data distribution and archival. However, these
data are also perfect for modeling the ionosphere using the mathematical
model outlined in the previous section.

*A. Data collection and cleaning*

For the results presented in this paper, data from day 193 (July 12) of
2002 were used, when CORS had 307 stations collecting data.

The data were collected from the CORS archive, in RINEX (Receiver
Independent Exchange) format. Then the raw carrier phases were
transformed into ∆TECS data, cleaned and sorted into tracks.

Each CORS station has one RINEX file for the 24 hour period of day 193,
year 2002. The data for each individual satellite were collected
together first, then points missing either L1 or L2 were thrown out. A
ten degree vertical angle cut off was used to maximize both crossovers
and spatial distribution and also so the impact of points low to the
horizon might eventually be studied. Tracks with fewer than 10
successive data epochs were considered too short to be useful and were
dropped.

For the remaining tracks, the ∆TECS values were computed using equation
9. These data then were investigated for cycle slips by discarding any
∆TECS values which were considered physically or statistically
unrealistic. For this study, "physically unrealistic" was set at 0.2
TECUs / km (Paul Spencer, NOAA, personal communication). Gradients
exceeding this value were fixed, if possible, using weighted averages of
neighboring gradients on the same track.

Finally, any loss of data for longer than 5 minutes constituted the end
of one track. The "cleaned" gradients, for the single CORS/Satellite
combination were then given unique track numbers, and the next
CORS/Satellite combination was similarly processed. All in all, the
number of "cleaned" tracks for this one day was 12,585, averaging about
1008 points (data epochs) per track.

[V. Initial Tests]{.smallcaps}

To test the rigor of the mathematical model of section III, small
tracknets (groups of tracks connected via crossovers) were formed and
their biases solved in a LSA. Then the modeled TECR values were
converted to double-difference (DD) values and compared to independently
computed double-differenced (DD) TECR values found through DD solutions
at the Geosciences Research Division of NGS (G. Mader, NOAA, personal
communication). A variety of LSA methods were tested. Presuming some
closed polygon could be found in the tracknet, then no a-priori
information was needed in the LSA, and the carrier-phase data could
entirely control the absolute level of TECR in the adjustment. However,
it was seen as useful to determine the impact of some a-priori
information, since this would be needed in any LSA of a tracknet where
no closed polygon of tracks could be formed, or even in a case where
a-priori information might be desired, though not necessary. This second
case presumes that the a-priori information is reliable and will add
useful information to the LSA.

A-priori information on the TECS biases for each track can be obtained
most quickly through a fitting of carrier phase TECS data to
pseudo-range TECS data. However, this presumes that the pseudo-range
data is available (not always true) and that the multi-path errors
average to zero (generally not true) so that the mean offset of phase
and pseudo-range would be an unbiased estimate of the true bias in the
phase-based TECR values. Considering the number of possible biases in
pseudo-range data, it is not totally surprising that the results using
only the carrier phase data generally ended up better than those that
introduced a-priori information from phase/pseudo-range fitting.

Figures 9 and 10 show two different examples of phase/pseudo-range
fitting. In Figure 9, one can clearly see that the phase data fits very
accurately through the highly noisy pseudo-range data; though,
considering how noisy the pseudo-range data are, it's difficult to know
what non-zero biases are hidden in that data set. In Figure 10, we can
see that the fit is even less clean and further calls into question the
usefulness of a phase/pseudo-range a-priori estimate.

![](media/image9.png){width="3.3743055555555554in"
height="2.501388888888889in"}

Fig. 9. A good fit between phase and pseudo-range based TECS values on
track 174.

![](media/image10.png){width="3.373611111111111in"
height="2.540277777777778in"}

Fig. 10. A poor fit between phase and pseudo-range based TECS values on
track 8366.

Many small (10-100 tracks) tracknets were formed and tested before
moving on to a larger adjustment. These tests served mostly to prove the
methodology, rather than provide significant insight into the ionosphere
itself. All results were similar, therefore only one example will be
presented here.

A pair of stations, viewing a pair of satellites were chosen to form a
double difference. In this example, they are CORS stations GODE and
RED1, receiving data from GPS SV1 and SV2. The four tracks (designated
4300, 4303, 9484 and 9487) formed by combining the 2 stations and the 2
satellites will be called the "DD tracks". It is the intention of this
small adjustment to investigate the best way to arrive at TECS values on
the 4 DD tracks. First off, it was necessary to find a small "net" of
tracks where they are all interconnected and which contains the 4 DD
tracks. Options in the adjustment are whether or not to have a 'closed
polygon' in the net, and whether to use none, some or all of the
a-priori bias estimates on the tracks in the net. Table 1 shows the
track numbers used in the tracknet.

[TABLE 1]{.smallcaps}

[Basic statistics of one small tracknet and 4 different]{.smallcaps}

[solutions for computing the Double Difference]{.smallcaps}

[TECR of tracks 4300/4303/9484/9487]{.smallcaps}

<table>
<colgroup>
<col style="width: 12%" />
<col style="width: 13%" />
<col style="width: 14%" />
<col style="width: 12%" />
<col style="width: 5%" />
<col style="width: 5%" />
<col style="width: 5%" />
<col style="width: 5%" />
<col style="width: 5%" />
<col style="width: 5%" />
<col style="width: 5%" />
<col style="width: 5%" />
</colgroup>
<tbody>
<tr class="odd">
<td></td>
<td></td>
<td></td>
<td></td>
<td colspan="4">Track used in solution?</td>
<td colspan="4"><p>A-priori</p>
<p>info used?</p></td>
</tr>
<tr class="even">
<td>Track #</td>
<td>CORS name</td>
<td>GPS SV #</td>
<td>Hours of day</td>
<td>1</td>
<td>2</td>
<td>3</td>
<td>4</td>
<td>1</td>
<td>2</td>
<td>3</td>
<td>4</td>
</tr>
<tr class="odd">
<td>4300</td>
<td>GODE</td>
<td>1</td>
<td>0.00 - 5.18</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td></td>
</tr>
<tr class="even">
<td>4303</td>
<td>GODE</td>
<td>2</td>
<td>0.00 - 5.93</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td></td>
</tr>
<tr class="odd">
<td>9484</td>
<td>RED1</td>
<td>1</td>
<td>0.04 - 5.15</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td>x</td>
</tr>
<tr class="even">
<td>9487</td>
<td>RED1</td>
<td>2</td>
<td>0.04 - 6.05</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td></td>
</tr>
<tr class="odd">
<td>2253</td>
<td>CONO</td>
<td>22</td>
<td>0.00 - 4.62</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td></td>
<td></td>
<td></td>
</tr>
<tr class="even">
<td>10146</td>
<td>SHK1</td>
<td>1</td>
<td>0.04 - 5.08</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td></td>
</tr>
<tr class="odd">
<td>11416</td>
<td>UIUC</td>
<td>25</td>
<td>0.00 - 2.70</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td></td>
<td></td>
<td></td>
</tr>
<tr class="even">
<td>12565</td>
<td>YOU1</td>
<td>20</td>
<td>0.04 - 2.42</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td></td>
</tr>
<tr class="odd">
<td>2224</td>
<td>COLB</td>
<td>25</td>
<td>0.00 - 2.91</td>
<td></td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td></td>
<td></td>
<td></td>
</tr>
<tr class="even">
<td>11580</td>
<td>UPO1</td>
<td>22</td>
<td>0.00 - 4.85</td>
<td></td>
<td>x</td>
<td>x</td>
<td>x</td>
<td></td>
<td>x</td>
<td></td>
<td></td>
</tr>
</tbody>
</table>

Four different solutions were investigated based on various options
mentioned above, but always with the goal being TECR values on the four
DD tracks. Note that the double difference exists from GPS time 0.04h to
5.15h, allowing over 5 hours of continuous data for testing. The first 8
tracks were chosen because they are the smallest set of tracks that all
connect one to another, such that the tracknet contains all 4 DD tracks.
However, these 8 tracks do not form a closed polygon on the ionosphere
shell; they contain only 7 crossovers and therefore would require
a-priori information to solve. The additional 2 tracks (2224, 11580)
were chosen because they combine with the first 8 tracks to form the
smallest tracknet that contains all 4 DD tracks and contains a closed
polygon. In fact, these 10 tracks contain 2 closed polygons[^1]. The 10
track tracknet contains 11 crossovers, and therefore could be solved
without any a-priori information.

Solution 1 contained the first 8 tracks (thus no polygons) and a-priori
estimates on all tracks that had pseudo-range data (6 of 8). Solution 2
contained all 10 tracks (thus 2 polygons) and a-priori estimates on all
tracks that had pseudo-range data (7 of 10). Solution 3 contained all 10
tracks and no a-priori estimates. Solution 4 contained all 10 tracks and
only one a-priori estimate (track 9484). These small tracknets were
adjusted and biases solved for. The formal (a-posteriori) error
estimates of the biases for each track, converted into cycles of delay
on L1 are presented in Table 2.

[TABLE 2]{.smallcaps}

[Statistics of the formal a-posteriori errors estimates of
track-by-track biases for the 4 solutions of DD 4300/4303/9484/9487.
Units are TECUs. (See also Table 1)]{.smallcaps}

  -------------- -------------- -------------- -------------- --------------
  Track \#       Solution 1     Solution 2     Solution 3     Solution 4

  4300           ± 3.5          ± 2.9          ± 0.1          ± 1.2

  4303           ± 8.8          ± 4.7          ± 0.2          ± 2.1

  9484           ± 9.3          ± 4.6          ± 0.2          ± 2.0

  9487           ± 9.4          ± 3.1          ± 0.1          ± 1.3

  2253           ± 13.6         ± 5.9          ± 0.3          ± 2.5

  10146          ± 9.7          ± 3.3          ± 0.1          ± 1.4

  11416          ± 6.5          ± 4.9          ± 0.2          ± 2.0

  12565          ± 6.1          ± 3.9          ± 0.2          ± 1.6

  2224           \-             ± 4.3          ± 0.2          ± 1.7

  11580          \-             ± 3.0          ± 0.1          ± 1.2
  -------------- -------------- -------------- -------------- --------------

One may draw some immediate conclusions from Table 2. The improvement
from solution 1 to solution 2 indicates that simply having a closed
polygon in the tracknet has stabilized the solution significantly.
Additionally, the improvement from solution 2 to solution 3 shows that
the a-priori information (based on pseudo-range/phase fitting) is
actually degrading the adjustment, causing dispersion in the adjusted
tracknet. Leaving the a-priori data out entirely in solution 3 seems to
strengthen the tracknet and yield a very rigid set of adjusted biases.
Finally, the slight degradation from solution 3 to solution 4 is due
solely to adding just 1 piece of a-priori information. This further
exemplifies the degrading nature of unnecessarily adding pseudo-range
data.

After these adjustments were performed, TECR values were converted to
cycles of delay on L1 for the entire track, and double differenced.
These were then provided to GRD for comparison to delay estimates that
came from an independent double-difference adjustment for the
coordinates of GODE and RED1 (Mader, ibid). The results are shown in
Figure 11.

![](media/image11.png){width="3.0625in" height="2.5618055555555554in"}

Fig. 11. Comparison between 4 different solutions for the DD TECR values
and a 5th independent solution from NGS positioning software showing
agreement between solutions 3 and 4 (yellow and dark blue) with the
independently computed DD TECR values (dark red).

As seen in Figure 11, solutions 3 and 4 have the best agreement to the
independent DD solutions, which further strengthens the argument that no
a-priori information needs to be put into the modeling of the
ionosphere. The agreement between solution 3 and the 5th (independent)
DD L1 delay values is at the 0.01 cycles value (RMS) and is remarkable,
considering TECR values were first computed as absolute values, and then
converted to DD values and absolutely no a-priori information was used
in solution 3. Numerous other examples, not detailed here, were
performed with similar conclusions.

[VI. Full Day Adjustment]{.smallcaps}

With the confidence that the small tracknet adjustment tests provided,
the next step was to attempt to compute a 24 hour adjustment of all
biases for all tracks that had at least one connection to a single,
master tracknet. Logical search operations were coded to find such a
tracknet. As it turns out, of the 12,585 tracks that were formed by CORS
on day 193 of 2002, there were 8298 which formed a master tracknet (that
is any of those 8298 tracks could be traced, via crossovers, to any
other of the 8298 tracks). The criteria for a crossover was that some
point on one track lay within a small space-time window \[0.1° x 0.1° x
60 seconds\] as a point on another track. Of those tracks that did not
connect to the master tracknet, most were either very short, very low on
the horizon, provided by a distant CORS station (e.g. Hawaii, Alaska,
Caribbean) or some combination of these factors.

Of the 8298 tracks (unknowns) that formed that master tracknet, there
were 16,896 crossovers (constraints) , thus providing a very redundant
system of linear equations. Sparse matrices were the obvious choice for
this adjustment. Each row of the A matrix corresponded to one constraint
(cross-over). The only non-zero elements in any given row of the A
matrix are in the 2 columns which correspond to the 2 tracks which make
the crossover. Although no immediate pattern of non-zero elements arises
from this situation, we nonetheless have only 0.024% of elements of A
that are non-zero.

Using the sparse matrix routines of NGS' "HEART OF GOLD" software \[5\],
this system of equations was solved for all 8298 track biases in 30
seconds, and the sparse profile elements of the dispersion matrix of the
adjusted biases following in a further 10 minutes of computations.

Once the adjusted biases and their dispersions were computed, we checked
the post-adjusted crossovers to see how well the master tracknet closed.
Table 3 shows the statistics of the post-fit crossovers. Table 4 shows
how well the formal error estimates of the adjusted biases came out.

[TABLE 3]{.smallcaps}

[Post-Fit Crossover statistics for master tracknet (TECUs)]{.smallcaps}

  ----------- ----------- ----------- ----------- ----------- -----------
  \#          Ave         STD         RMS         Min         Max

  16,896      -0.004      ±0.51       0.51        -3.7        +4.0
  ----------- ----------- ----------- ----------- ----------- -----------

[TABLE 4]{.smallcaps}

[Statistics of A-Posteriori Standard Deviations for post-fit master
tracknet biases (TECUs)]{.smallcaps}

  ----------------- ----------------- ----------------- -----------------
  \#                Ave STD           Min STD           Max STD

  8298              ±1.1              ±0.22             ±10.7
  ----------------- ----------------- ----------------- -----------------

In addition to these statistics, it is useful to consider the actual
TECR values for every point on every track in the master tracknet.
Applying the adjusted TECS bias of each track to the known ∆TECS values,
and then generating TECR values, the following statistics are found in
Table 5.

[TABLE 5]{.smallcaps}

[Post-Fit TECR values for entire 24 hour master tracknet
(TECUs)]{.smallcaps}

  ------------- ----------- ----------- ----------- ----------- -----------
  \#            Ave         STD         RMS         Min         Max

  10,613,013    14.90       ±6.03       16.07       -7.5        +58.4
  ------------- ----------- ----------- ----------- ----------- -----------

As seen in Table 5, the adjusted TECR values are reasonable (on average)
but with poor outliers (below zero or above 50 TECUs). The negative TECR
values, however, are rare, (0.08% of all TECR values) as are those above
50. Aside from these outliers, we can show the generally believable
nature of the ionosphere more clearly. The adjusted TECR values were
averaged in 2 hour intervals and plotted over 24 hours. Additionally,
values for TECR over the same geographic area from 5 IGS (International
GPS Service) computational centers (COD, EMR, ESA, JPL and UPC) were
downloaded from their archive and plotted. See Figure 12.

![](media/image12.png){width="3.372916666666667in"
height="2.6305555555555555in"}

Fig. 12. The 5 IGS models compared to the method of this paper showing
good agreement at night and early morning

Note that the agreement is much better between this method and the
global models at night and morning (central parts of plot). The
disagreements for the late afternoon (left side of plot) and late
morning (right side of plot) could be due to a number of sources. First,
this method used over 300 stations just in the USA, while all 5 global
models averaged only 120 stations for the entire globe. Additionally,
the geographic and temporal spacing for the global models was 5° x 2.5°
x 2 hours, whereas this method models the ionosphere at every data epoch
(1, 5, 15 or 30 seconds) on the densely packed tracks themselves.
Additionally, the global models assume an ionosphere shell height of 450
km while this study used 300 km. Because the cut-off angle used in this
paper was 10 degrees, the errors in mapping TECS to TECR via equation 1
can lead to errors in remote tracks, such as those over the ocean areas
where only a few (shoreline) stations can see the rising satellites
through the ionosphere, and then only low on the horizon. Finally, the
disagreement occurs during the most active time of the ionosphere itself
as well, which also propagate into the errors of using equation 1.

Considering the good agreement between this method and the 5 global
methods during at least the night/morning half of the day, and the
explainable (and possibly repairable) sources of disagreement during the
rest of the day, Figure 12 is encouraging as a first step. This is
especially so when one considers that no pseudo-range fitting or other
a-priori methods of solving for the biases were used. At the very least,
the method is well-suited to double difference ionosphere modeling, and
shows great promise as a method for absolute ionosphere modeling,
provided further study can be performed on the disagreements of Figure
12.

At this point, the ionosphere has been modeled, but at discrete times
(data epochs) and discrete locations (tracks). We then turned our
attention to the question of how best to use such data, and
subsequently, how best to distribute it. Clearly, visualization of the
ionosphere would be an aid to anyone interested in understanding or
using this model. That would generally imply transforming our values
onto a grid. However, one of the ultimate goals of this project was to
provide, at the best possible accuracy, values of ionospheric delay to
surveyors using a GPS receiver in the field. Going from tracks to a grid
back to a surveyor's track would mean two sources of interpolation
error. Going directly from the adjusted master tracknet to a surveyor's
track would seem to be a more accurate path.

Ultimately, both were investigated, and their advantages and
disadvantages are discussed in the next sections.

[VII. Gridding]{.smallcaps}

Putting aside for the moment the singular application of providing
ionosphere delays to individual surveyors, it would seem that the
obvious path to provide our ionosphere model to the general public would
be in the form of a grid of TECR values at discrete points in time. Such
a grid lends itself to easy distribution, interpolation, analysis,
conversion to other units, comparison to other ionosphere models and of
course, visualization. The problem of transforming random points to a
grid has a long history. Splines, collocation and polynomials all have
advantages and disadvantages. What was seen as important here was
finding a method that would allow the spatial and temporal connections
to remain in the grids, while allowing fast and accurate gridding. The
method chosen for this paper was one of many possible methods and is
certainly not meant to be the final choice. However, for the sake of
expediency, a few quick experiments lead to the following scheme:

Step 1: Fit a 3-D polynomial of degree and order 1 (called TECR0) to the
TECR data:

$\text{TECR}0(\varphi,\lambda,t) = \sum_{i = 0}^{1}{}\sum_{j = 0}^{1}{}\sum_{k = 0}^{1}{a_{i,j,k}\varphi^{i}}\lambda^{j}t^{k}$
(21)

Step 2: Remove TECR0 from TECR

Step 3: Block-mean the residual data in blocks of size 0.1° x 0.1° x 15
min

Step 4: Transform block-means to 2-D gridded values every 15 minutes
using splines in tension \[6\], and a tension parameter of T = 0.1

Step 5: Restore TECR0 to gridded residuals

and, for the sake of visualization only, to avoid serious edge effects:

Step 6: Mask the final grids by using a 100km buffer around all φ / λ
locations which had at least one pierce point over the 24 hour period.

One could chose many different versions of this scheme (different
polynomial, or using an a-priori grids (such as the IGS global grids); a
different tension parameter; a different size for block-means) or a
different scheme altogether (e.g. Fourier series; Collocation). The
emphasis of this study is not to study which of these schemes is best --
that will come in a future paper. The emphasis is to show that even the
simplest scheme can yield results that are both visually pleasing and
statistically reasonable.

Figure 13 shows the final, gridded and masked ionosphere values for year
2002, day 193, at time 13500 sec. This is one of 97 grids which make up
the entire 24 hour period.

![](media/image13.png){width="3.3229166666666665in" height="2.25in"}

Fig. 13. Modeled TECR data (masked to remove edge effects) at t=13500
sec, with raw tracks shown for the 15 minute period surrounding this
time.

All 97 grids were animated, and made available at
<http://www.ngs.noaa.gov/IONO/f.gif>. Such animation serves two
purposes. It is an invaluable tool for analysis of actual ionospheric
activity, but also allows a very quick method of determining outliers.
We are currently in the process of re-cleaning the data based on this
animation and the information it yielded. It should be noted that the
resolution, in both space and time, of this ionosphere grid is much
better than anything previously seen in global models.

Once the grids were computed, the question arose of how well we could go
from a grid back to a track. (This would obviously be the way our grids
would be used). Again, we adopted the strategy of testing the simplest
solution first, knowing that more elaborate schemes of interpolation
should only improve upon the initial results. Therefore, taking the
grids, and using a basic tri-linear interpolation of the closest 8
points, we attempted to re-create the original adjusted master tracknet.
The statistics of that experiment are shown in Table 6.

[TABLE 6]{.smallcaps}

[Misfit of original tracknet to tracknet re-generated from gridded data
(TECUs)]{.smallcaps}

  --------------- ------------- ------------- ------------- -------------
  \#              Ave           STD           Min           Max

  10,613,013      +0.001        ±0.38         -7.8          +8.8
  --------------- ------------- ------------- ------------- -------------

As can be seen in Table 6, the ability to go from tracknet to grid and
back to tracknet can be done, in a very simple method to almost perfect
accuracy and with a precision of 0.38 TECUs. Certainly such a precision
could be improved with better gridding and interpolation algorithms, but
on the other hand, even the best of these methods may not improve the
ability to correctly predict the ionosphere in locations removed in
space and time from the data of the master tracknet. Such areas arise
when GPS receivers are outside of (but near) the USA, or even in the USA
itself in areas where CORS is sparse (such as the north central region;
see Figure 8).

Given that the gridded data appears to be a useful for visualization and
even a reasonable method of accessing the information contained on the
tracks themselves, we turn our attention next to the possibility of
accurately predicting ionosphere data on a track by directly
interpolating off of the master tracknet itself.

[VIII. Interpolation Via a TIN]{.smallcaps}

It was already seen that, in addition to being a good visual tool, a
gridded ionosphere model can be a reasonable predictor, at least to a
precision of 0.38 TECUs. It was hoped that by skipping the gridding
altogether, an even greater precision could be achieved. The scheme
relies on a method quite like that of a Triangulated Irregular Network
(TIN).

Consider 3 CORS stations that form a triangle; inside this triangle is a
4th GPS receiver. If, at the same epoch all four of these receivers are
receiving data from the same GPS satellite, and if the 3 CORS stations
can compute the ionosphere through the adjustment of the master
tracknet, it stands to reason that we should be able to estimate what
the ionosphere seen by the 4th GPS receiver should be, based on the 3
surrounding stations.

To test this hypothesis, each CORS station contributing to the master
tracknet was investigated as follows: First, a variety of "surrounding
triangles" (if any) were found, based on all other CORS stations in the
master tracknet. Four types of triangles were studied to see if the
choice of triangle would have a noticeable advantage over the other
types. The four types of surrounding triangles are outlined in Table 7.

[TABLE 7]{.smallcaps}

[Different triangle types for TIN experiments]{.smallcaps}

  ---------------- ------------------------------------------------------
  Triangle Type    Description

  1                with the smallest area

  2                with the smallest maximum side

  3                with the smallest maximum distance from CORS to
                   central receiver

  4                with the largest minimum angle
  ---------------- ------------------------------------------------------

When the triangles were defined, the adjusted master tracknet was run
through the comparison program. Basically it would make sure that all 4
CORS receivers were seeing the same satellite at the same epoch. It
would then use the TECR values from the 3 corner CORS stations to
approximate the TECR value at the central CORS station, using the
geometry of the triangle on the ground. (An alternate, and possibly more
robust, solution would be to re-compute the triangular geometry formed
by the 4 pierce points at every epoch, but this has the disadvantage of
needing the relations between all 4 stations re-computed every epoch).
This interpolated TECR value was then compared with the actual value of
TECR from the central CORS site as seen in the master tracknet. The
results are shown in Table 8.

[TABLE 8]{.smallcaps}

[Statistics of track-regeneration via TIN methods, based on 4 triangle
types (TECUs)]{.smallcaps}

  ------------ -------------- ----------- ----------- ---------- -----------
  Triangle     \# points      Ave         STD         Min        Max
  Type                                                           

  1            1,950,266      -0.010      ±0.28       -10.8      +11.8

  2            2,761,325      -0.005      ±0.25       -11.1      +11.6

  3            2,569,277      -0.006      ±0.25       -11.1      +12.0

  4            1,933,781      -0.029      ±0.29       -10.2      +11.6
  ------------ -------------- ----------- ----------- ---------- -----------

A few important things can be gleaned from Table 8. First, not every
CORS station in the master tracknet could be surrounded by 3 other such
CORS stations. As such, those data were not part of this analysis.
Additionally, since CORS stations can have data rates of 1, 5, 15 or 30
seconds, it is possible for significant data to be skipped in this
analysis if one of the CORS stations has a lower data rate than the
other 3, since common epochs were a requirement for comparison.

Finally, it should be noted that independent of the type of triangle
chosen, this method yields results that are around 0.25 TECUs, which is
an improvement over the gridding method of section VIII.

One must be careful not to presume that the ability of either the
gridding or TIN methods to re-create the original TECR data is a measure
of absolute accuracy. Certainly both the gridding and TIN methods have
high precision, but the accuracy of the method itself will only come
through independent verification of this model. Currently, the best
comparisons available are from coarse (in both space and time) global
models, but the difficulties with such a comparison have already been
outlined.

In the future, actual application of this method toward directly
reducing the time needed to resolve integer ambiguities will be the best
test and guidance for future directions.

[IX. Conclusions and Future Directions]{.smallcaps}

This paper was an introduction to a new method for unambiguously
computing the total electron content of the ionosphere using ambiguous
carrier phase GPS data only. While certain assumptions were required for
this method to work, notably the use of the constant-height shell model,
it nonetheless was a mathematically sound method. There were no tricks,
iterations or other ways to get around the ambiguity of carrier phase
data. Biases were found during daylight hours between this model and
global models at the 5+ TECU level. But these appear to be solvable with
future research. Aside from these biases from the global model, the
method in this paper appears to have a precision approaching 0.25 TECUs.
By converting TECUs into cycles of delay on L1 and L2, we can say that
the method described in this paper has the ability to accurately model
the ionosphere to a precision of 0.21 cycles on L1 (\~ 4.0 cm) and 0.27
cycles on L2 (\~ 6.6 cm). When double-differenced, a precision of 0.01
TECU was found.

The adjustment of thousands of crossovers for the solution of thousands
of biases in only 30 seconds meant that this method can be easily
adapted to a sequential least-squares method, where old tracks can be
updated and new tracks allowed to enter the master tracknet and
sequential adjustments performed in sequential epochs. All that was
required for a track to yield valid ionosphere data was that a single
crossover must exist between a track and the master tracknet. Once the
single bias for the track is known, all TECR values along that track are
therefore determined in absolute value. Whether this ionosphere model is
used in a

grid or a TIN method, its speed of computation, precision, and
dependence on an existing GPS network make it a powerful new tool for
ionospheric information distribution. But serious questions remain, all
of which are slated to be covered in future studies. Questions to be
answered include: How many CORS stations are needed before no further
ionospheric information is gained? What is the best method to test and
validate the data? How does this model behave on days of significant
ionospheric activity? How can we resolve the existing discrepancies with
the global models? And, of greatest interest to geodesists, what
improvement is gained in integer ambiguity resolution?

These questions are all valid, but the primary question of this study
has been answered: It is mathematically possible, almost trivial, to
compute the unambiguous value of TEC in the ionosphere from ambiguous
carrier phase GPS data only. The next phase is testing this method for
applicability to the integer ambiguity resolution problem, and
discovering why biases exist between this method and the global models.

[Acknowledgment]{.smallcaps}

The author wishes to acknowledge those who have provided input to this
process along the way. To Dr. Gerry Mader for providing the initial
encouragement to seek out a way to model the ionosphere from the CORS
network itself, and to Dr. Gerry Mader and Dr. Dennis Milbert for a
continuing series of technical discussions throughout this entire
research project.

[References]{.smallcaps}

\[1\] Odjik, D., 2000: Weighting ionospheric corrections to improve fast
GPS positioning over medium distances, Proceedings of the Institute of
Navigation 2000 meeting

\[2\] Teunissen, P.J.G. and D. Odjik, 2003: Rank-defect integer
estimation and phase-only modernized GPS ambiguity resolution, Journal
of Geodesy, v. 76, pp. 523-535.

\[3\] Wanninger, L. and M. Manja, 2001: Carrier-phase multipath
calibration of GPS reference stations, NAVIGATION, v. 48, n. 2, pp.
113-124.

\[4\] Langley, R.B., 1996: Propogation of the GPS Signals, Lecture Notes
in Earth Sciences, v. 60 (Kleusberg and Teunissen Eds.)

\[5\] Milbert, D.G.,1984: HEART OF GOLD: Computer routines for large,
sparse, lease squares computations. NOAA Technical Memorandum NOS
NGS-39, 1984.

\[6\] Smith, W. and P. Wessel, 1990: Gridding with continuous curvature
splines in tension, Geophysics, v. 55, n. 3, pp. 293-305.

[^1]: The two polygons are 8 sided and 3 sided, respectively, and have
    the following tracks as sides:
    9484-2253-10146-11580-2224-4300-11416-12565 and 9487-10146-11580.