GNSS satellite orbit determination and station network analysis
This cookbook chapter describes an example of global GNSS processing as done by analysis centers of the International GNSS Service (IGS). Resulting products usually comprise:
- Satellite orbits, clocks, and signal biases
- Station positions, clocks, signal biases, and troposphere estimates
- Earth orientation parameters
Scientific details about the underlying processing approach and the applied parametrizations, models, and corrections can be found in a doctoral thesis available under DOI 10.3217/978-3-85125-885-1.
An example scenario for this task is available at https://ftp.tugraz.at/outgoing/ITSG/groops/scenario/scenarioGnssNetwork.zip. It includes GROOPS scripts and data for the example, but not the general GROOPS data and metadata found at https://ftp.tugraz.at/outgoing/ITSG/groops (data folder or zipped archive). The scenario generally represents what is described in this cookbook, but may slightly differ in certain settings.
Note: Global GNSS processing can become very computationally intensive. Depending on the number of satellites and stations, the observation and processing sampling, and parametrizations it can quickly exceed the capabilities of a normal desktop computer and may require computer clusters or number crunchers (see section Parallelization).
Data preparation
Most of the required metadata files are provided in GROOPS file formats at https://ftp.tugraz.at/outgoing/ITSG/groops. These files are regularly updated.
Data that has to be gathered from other sources comprises:
- Receiver observations: GNSS measurements converted from RINEX format (see RinexObservation2GnssReceiver).
- Approximate orbits: broadcast or precise orbits in CRF for orbit integration (see GnssRinexNavigation2OrbitClock or Sp3Format2Orbit).
- Approximate clocks: broadcast or precise clocks (see GnssRinexNavigation2OrbitClock or GnssClockRinex2InstrumentClock)
The example scenario includes a small set of this data.
The script 010groopsConvert.xml
can be used to convert these external formats into GROOPS formats.
Prepare a station list file that contains the stations to be processed. Each line can contain more than one station. The first station in each line that has data available is used for the processing. If your network contains more than 60-70 stations, it is recommended to start processing with a core network (see Advanced). In this case, define an additional core station list file that can also have multiple stations per line.
Preprocessing: Orbit integration
Numerical integration of the satellite orbits is the first step in global GNSS processing. Dynamic orbits are integrated based on force models and then fitted to the approximate orbits by estimating their initial state and additional empirical parameters for solar radiation pressure to improve the orbit fit. The resulting variational equations file contains the integrated orbit, derivatives with respect to the satellite state vector, attitude, Earth rotation and satellite model.
Orbit preprocessing is covered by the script 020groopsGnssPreprocessing.xml
in the example scenario.
It is recommended to perform the steps below in a loop over all
satellites/PRNs using LoopPrograms. To get the relation between {prn}
and {svn}
setup
an additional loop:platformEquipment inside
loop:loop with
- inputfilePlatform: the old inputfileTransmitterInfo
-
equipmentType =
gnssAntenna
-
variableLoopName =
block
-
variableLoopSerial =
svn
-
variableLoopTimeStart =
svnTimeStart
-
variableLoopTimeEnd =
svnTimeEnd
-
condition:expression
-
expression =
(svnTimeStart <= loopTime) && (loopTime < svnTimeEnd)
{prn}
:
- InstrumentResample: resample approximate orbits from data preparation to target sampling (e.g., 1 minute) by defining a timeSeries based on a method:polynomial (polynomialDegree=
7
, maxDataPointRange=7200
, maxExtrapolationDistance=900
). -
OrbitAddVelocityAndAcceleration: add velocity via running polynomial (polynomialDegree=
2
) derivation (needed for attitude computation) - SimulateStarCameraGnss
-
PreprocessingVariationalEquation:
-
inputfileSatelliteModel=
{groopsDataDir}/gnss/transmitter/satelliteModel/satelliteModel_boxWing.{svn}.xml
- inputfileOrbit: the resampled approximate orbit from InstrumentResample
- inputfileStarCamera: the attitude file from SimulateStarCameraGnss
- forces: see below
-
gradientfield:potentialCoefficients: a static gravity field (e.g. GOCO06s) with maxDegree=
4
.
Force models usually include:
- gravityfield:potentialCoefficients: static gravity field (e.g. GOCO06s)
-
gravityfield:trend
- gravityfield:potentialCoefficients: trend component of time-variable gravity field (e.g. GOCO06s)
60
is more than enough.The result of the preprocessing should be a variational equations file, a reduced dynamic orbit file from PreprocessingVariationalEquationOrbitFit and an attitude file from SimulateStarCameraGnss for each satellite.
GNSS processing
The script 030groopsGnssProcessing.xml
in the example scenario
implements the following steps and settings.
These are the settings for GnssProcessing. If not otherwise stated use the default values.
The first step is setting the processing sampling, in this example it is 30 seconds.
The processing interval usually is a single 24-hour day, therefore define
timeSeries:uniformSampling with timeStart=<mjd>
,
timeEnd=<mjd>+1
, sampling=30/86400
(processing sampling).
Add the appropriate transmitters:gnss (e.g. GPS, GLONASS, and Galileo) and provide the required files:
- inputfileOrbit from preprocessing
- inputfileAttitude from preprocessing
- inputfileClock from data preparation
The following settings are needed in receiver:stationNetwork:
- inputfileStationList: list of all stations to be processed
- inputfileObservations: The converted RINEX observation files.
- tidalDisplacement: Use the settings described in receiver:stationNetwork.
-
excludeType: Signals you might want to exclude are
C*?G
(old unknown GPS code observations),*3*R
(GLONASS G3 freq.),*6*E
(Galileo E6 freq.).
Add the following parametrizations and define the outputfiles you are interested in inside each of them:
- ionosphereSTEC: add a constraint of sigmaSTEC=
40
-
ionosphereMap: add temporal:splines
with linear (degree=
1
) 2-hourly splines -
clocks: optionally change selectTransmitterZeroMeans to wildcard with name=
G*
to align clocks to mean over GPS (instead of all) satellites - signalBiases
- ambiguities
- codeBiases
- tecBiases
-
temporalBias: time-variable GPS L5 phase bias with type=
L5*G
and parametrizationTemporal:splines with degree 3 and hourly nodes. - staticPositions
-
troposphere: select troposphere:viennaMapping
with the appropriate vmf3grid file. Add troposphereWetEstimation:splines
with linear (degree=
1
) 2-hourly splines and troposphereGradientEstimation:splines with linear daily splines. -
transmitterDynamicOrbit:
provide inputfileVariational=
preprocessing/variational.{prn}.dat
from the preprocessing step. Add parametrizationAcceleration:gnssSolarRadiation and stochasticPulse:irregular parameter at center of day to further improve orbit fit. -
earthRotation
- estimatePole:constant: polar motion
-
estimatePole:trend polar motion rate
(at center of day with timeStep=
1
) -
estimateUT1:trend: length of day (at center of day with
timeStep=
-1
to match IGS sign convention)
signalBias.L5*
, sigma=5
meters, and relativeToApriori=yes
troposphere*
, sigma=5
meters, and relativeToApriori=yes
stochasticPulse*
, sigma=0.1
micrometers/second.
Finally, define the processingSteps. This can be overwhelming at first, but offers a lot of flexibility. The example script uses a 5-minute processing sampling with subsequent clock densification to 30 seconds.
- selectEpochs: with nthEpoch=
10
to reduce sampling to 5 minutes. -
selectParametrizations:
disable
constraint.STEC
,*VTEC
,*.tecBiases
as the ionosphere parameters are estimated in the final steps only. -
estimate: with maxIterationCount=
6
- resolveAmbiguities
-
selectParametrizations:
enable
*
(all) parameters -
estimate: with maxIterationCount=
4
: final iterations (with 5-minute sampling) and ionosphere parameters -
group: clock densification to 30-second sampling
-
selectEpochs: with nthEpoch=
1
to set full 30-second sampling -
selectParametrizations:
disable
*
(all) parameters and reenable*.clock*
and*.STEC
parameters -
estimate: with maxIterationCount=
6
-
writeResults: with suffix=
30s
to write 30-second clock files
With some additional steps, the full 30-second sampling can be used to estimate all parameters (not only the clocks). These steps are disabled in the example script, as they require at least 16 GB of system memory. In this case, it is not necessary to separately write the 30-second clock files as listed above.
- selectEpochs: with nthEpoch=
1
to set full 30-second sampling -
selectNormalsBlockStructure: As the system of normal
equations can be very large, the memory consumption might be reduced with keepEpochNormalsinMemory=
no
. In this case the epoch parameters are directly eliminated during the accumulation and reconstructed in the solving step. This might lead to longer computation times. -
estimate: with maxIterationCount=
2
: final iterations with full sampling
Advanced: Processing large station networks
Processing large station networks requires some additional steps to keep the computational load to a reasonable degree. The general processing strategy is to first process a well-distributed subset of stations (i.e. a core network) to get good estimates of all satellite parameters, which then enables integer ambiguity resolution (IAR). Once the ambiguities of the core network are resolved and stable estimates for satellite phase biases are available, all other (non-core) stations can be processed individually (including IAR) while keeping the satellite parameters fixed. At last, all stations can be processed together with all satellite parameters and ionosphere parameters.
Let's start with the processingSteps of the core network:
- selectReceivers with selectReceivers:file using the core network station list file from data preparation as inputfileStringTable.
-
selectEpochs: with nthEpoch=
10
to reduce sampling to 5 minutes. -
selectParametrizations:
disable
constraint.STEC
,*VTEC
,*.tecBiases
as the ionosphere parameters are estimated in the final steps only. -
estimate: with maxIterationCount=
6
- resolveAmbiguities
-
estimate: with maxIterationCount=
4
: final iterations (with 5-minute sampling)
Now all other (non-core) stations can be processed separately:
- forEachReceiverSeparately:
with selectReceivers:file inside selectReceivers:exclude
using the station list from the core network above to process all non-core stations individually with fixed transmitter parameters
-
processingStep:estimate: with maxIterationCount=
6
- processingStep:resolveAmbiguities
-
processingStep:estimate: with maxIterationCount=
4
Next all stations are processed together with all parameters:
- selectReceivers: with selectReceivers:all
-
selectEpochs: with nthEpoch=
1
to set full 30-second sampling -
group: clock densification to 30-second sampling
-
selectParametrizations:
disable
*
(all) parameters and reenable*.clock*
and*.STEC
parameters -
estimate: with maxIterationCount=
6
*
(all) parameters.
no
4
: final iterations with full sampling and all parameters