GNSS Measurement Model¶

Purpose¶

This specification defines the GNSS receiver measurement model implemented by GnssMeasurement and GNSSMeasurements. The model generates pseudorange, Doppler, and carrier phase observations from receiver state, GNSS transmitter ephemerides, receiver/transmitter clock corrections, and optional propagation delays.

State, Time, and Frame Contract¶

All vectors in this specification are expressed in options.frame. The current C++ implementation requires receiver and transmitter states to share a defined common inertial frame; the lunar filtering project uses Frame::MOON_CI for receiver states and GNSS channel states after constellation setup.

The receiver state vector is interpreted through GnssMeasurementOptions::indices:

\[x = \begin{bmatrix} \cdots & r_R^\mathsf{T} & \cdots & v_R^\mathsf{T} & \cdots & b_R & d_R & \cdots \end{bmatrix}^\mathsf{T},\]

where

\[r_R \in \mathbb{R}^3 \quad [\mathrm{m}], \qquad v_R \in \mathbb{R}^3 \quad [\mathrm{m/s}].\]

The receiver clock bias b_R and drift d_R are represented in options.clock_bias_unit. They are converted to SI seconds before entering the range equations:

\[\Delta t_R = \operatorname{seconds}(b_R), \qquad \dot{\Delta t}_R = \operatorname{seconds}(d_R).\]

Epochs passed to Compute and Precompute are receiver signal-reception epochs:

\[t_R \equiv t_\mathrm{receive}.\]

The scale of t_R is options.receive_time_scale. GNSS constellation ephemerides are tabulated in TAI seconds and are evaluated in options.ephemeris_time_scale; the current implementation enforces TAI for the ephemeris scale. Cislunar simulations commonly set the receiver receive epoch scale to TDB and let the measurement model convert receive epochs to TAI for GNSS ephemeris interpolation.

GNSS Transmitter State Interpolation¶

For each PRN, GnssConstellation::SetSatelliteStates receives tabulated states

\[\{(t_k,\; r_S(t_k),\; v_S(t_k))\}_{k=0}^{N-1}.\]

The constellation fits a piecewise Chebyshev model for the six-dimensional state:

\[\begin{split}y_S(t) = \begin{bmatrix} r_S(t) \\ v_S(t) \end{bmatrix} \approx \sum_{n=0}^{p} c_n T_n(\tau),\end{split}\]

where

\[\tau = \frac{t - t_\mathrm{mid}}{t_\mathrm{radius}}, \qquad -1 \le \tau \le 1.\]

Runtime satellite state evaluation uses this Chebyshev model. The original state table may still be carried in each GnssChannel as ephemeris-message metadata.

Light-Time Iteration¶

For a receiver position r_R(t_R) and transmitter PRN s, the transmit epoch is solved iteratively in the ephemeris time scale:

\[t_S^{(0)} = t_R,\]

\[t_S^{(i+1)} = t_R - \frac{ \left\lVert r_R(t_R) - r_S(t_S^{(i)}) \right\rVert }{c}.\]

The iteration stops when

\[\left| t_S^{(i+1)} - t_S^{(i)} \right| < \epsilon_\mathrm{lt},\]

or when options.light_time_max_iterations is reached. The convergence tolerance is options.light_time_tolerance_s.

Clock terms and group delays are not used to solve the physical light-time epoch; they enter the measurement equations below.

Geometry¶

After light-time convergence, define

\[r_S = r_S(t_S), \qquad v_S = v_S(t_S),\]

\[\Delta r = r_R - r_S, \qquad \Delta v = v_R - v_S,\]

\[\rho = \lVert \Delta r \rVert, \qquad \hat{u} = \frac{\Delta r}{\rho},\]

\[\dot{\rho} = \frac{\Delta r^\mathsf{T}\Delta v}{\rho}.\]

Transmitter Clock Corrections¶

The transmitter effective clock bias is

\[\Delta t_S^\mathrm{eff} = \Delta t_S + \Delta t_S^\mathrm{rel} - \Delta t_S^\mathrm{grp},\]

where tx_clock_bias_s is \(\Delta t_S\), relativistic_correction_s is \(\Delta t_S^\mathrm{rel}\), and group_delay_s is \(\Delta t_S^\mathrm{grp}\).

The current transmitter relativistic correction is the state-based broadcast GNSS form

\[\Delta t_S^\mathrm{rel} = -\frac{2 r_S^\mathsf{T} v_S}{c^2}.\]

The transmitter effective clock drift is

\[\dot{\Delta t}_S^\mathrm{eff} = \dot{\Delta t}_S.\]

Propagation Delays¶

Shapiro Delay¶

When options.apply_shapiro_delay is true, the one-way gravitational delay from the Sun is represented as a range correction:

\[\Delta \rho_\mathrm{Shapiro} = \frac{2\mu_\odot}{c^2} \ln \left( \frac{ r_R^\odot + r_S^\odot + \rho }{ r_R^\odot + r_S^\odot - \rho } \right),\]

where

\[r_R^\odot = \lVert r_R - r_\odot(t_R) \rVert, \qquad r_S^\odot = \lVert r_S - r_\odot(t_R) \rVert.\]

The Sun position \(r_\odot(t_R)\) is supplied by SetSunPositionProvider or by the default Sun provider. The legacy options.shapiro_mu and options.shapiro_body_position fields are not used by the current implementation.

Ionosphere / Plasma Delay¶

The ionosphere/plasma correction is represented as a positive code delay in meters:

\[\Delta \rho_\mathrm{plasma} \ge 0.\]

It is disabled by default. When options.apply_ionosphere_plasma_delay is true, the delay is supplied by one of:

SetCustomIonospherePlasmaDelayModel for a custom online link-by-link model.
SetIonospherePlasmaRayTraceOptions for built-in link-by-link evaluation with the plasma ray tracer in environment/plasma/tec/raytrace.
SetBatchCustomIonospherePlasmaDelayModel for a custom precompute model over all epochs and all channels.
options.default_ionosphere_plasma_delay_m when no provider is installed.

The built-in raytrace option calls pecsim::trace_ray directly. The caller must supply the pecsim::RayTraceConfig and, when desired, the IRI model selection; GNSSMeasurements does not choose Kp, step size, integrator, correction method, cutoff radius, or GCPM/IRI settings. The GNSS channel frequency is copied into RayTraceConfig::freq_Hz by default, but this can be disabled through use_channel_frequency. Receiver and transmitter positions are converted from options.frame to GnssIonospherePlasmaRayTraceOptions::raytrace_frame before tracing, and the receive epoch is converted from options.receive_time_scale to raytrace_epoch_scale.

The raytracer returns a path profile. The GNSS measurement uses PathProfile::tec_delay_m by default:

\[\Delta \rho_\mathrm{plasma} = \frac{40.3}{f^2}\int_\Gamma n_e \, ds .\]

When delay_mode is TEC_PLUS_HIGHER_ORDER, the second- and third-order terms reported by the raytracer are added to this code delay.

The sign convention is dispersive:

\[\text{code range uses } +\Delta \rho_\mathrm{plasma}, \qquad \text{carrier phase range uses } -\Delta \rho_\mathrm{plasma}.\]

Observable Equations¶

Let

\[\lambda = \frac{c}{f}.\]

The receiver clock range and range-rate terms are

\[\Delta \rho_R = c \Delta t_R, \qquad \Delta \dot{\rho}_R = c \dot{\Delta t}_R.\]

The transmitter clock range and range-rate terms are

\[\Delta \rho_S = c \Delta t_S^\mathrm{eff}, \qquad \Delta \dot{\rho}_S = c \dot{\Delta t}_S^\mathrm{eff}.\]

Pseudorange¶

The pseudorange observable is

\[P = \rho + \Delta \rho_\mathrm{Shapiro} + \Delta \rho_\mathrm{plasma} + \Delta \rho_R - \Delta \rho_S.\]

Doppler¶

The Doppler observable in Hz is

\[D = - \frac{ \dot{\rho} + \Delta \dot{\rho}_R - \Delta \dot{\rho}_S }{\lambda}.\]

The current implementation does not include time derivatives of Shapiro or plasma delay in Doppler.

Carrier Phase¶

The non-integer carrier phase observable in cycles is

\[\Phi = \frac{ \rho + \Delta \rho_\mathrm{Shapiro} - \Delta \rho_\mathrm{plasma} + \Delta \rho_R - \Delta \rho_S }{\lambda} + \phi_0,\]

where \(\phi_0\) is phase_bias_cycles.

If the measurement state includes an integer ambiguity \(N\), the emitted carrier observable is

\[\Phi_\mathrm{meas} = \Phi + N.\]

Visibility and Channel Selection¶

Visibility is receiver-dependent and is evaluated in GNSSMeasurements. For each occluding body, the line of sight is retained only if the spherical body does not block the receiver-transmitter segment.

For near-surface points, visibility is evaluated using an elevation threshold

\[e > -10^\circ.\]

For space-to-space links, visibility is evaluated using a tangent-line occlusion test. C/N0 thresholding is then optionally applied:

\[C/N_0 > (C/N_0)_\mathrm{min}.\]

Noise and Covariance¶

For each selected channel, the measurement covariance is diagonal in the current implementation:

\[R = \operatorname{diag} \left( \sigma_P^2, \sigma_D^2, \sigma_\Phi^2 \right),\]

using the subset and order requested by options.observables.

When C/N0 is finite, the tracking-loop noise models provide:

\[\sigma_P \quad [\mathrm{m}], \qquad \sigma_D \quad [\mathrm{Hz}], \qquad \sigma_\Phi \quad [\mathrm{cycles}].\]

If a channel does not have a finite positive standard deviation for an observable, the filter wrapper uses a fallback standard deviation of 1 in that observable’s native unit.

Jacobian Convention¶

The analytic Jacobian currently includes first-order derivatives of geometric range, range rate, receiver clock bias, receiver clock drift, and optional carrier integer ambiguity with respect to the receiver state.

For pseudorange:

\[\frac{\partial P}{\partial r_R} = \hat{u}^\mathsf{T}, \qquad \frac{\partial P}{\partial b_R} = c \frac{\partial \Delta t_R}{\partial b_R}.\]

For Doppler:

\[\frac{\partial D}{\partial r_R} = -\frac{1}{\lambda} \left( \frac{\Delta v}{\rho} - \frac{\dot{\rho}\Delta r}{\rho^2} \right)^\mathsf{T},\]

\[\frac{\partial D}{\partial v_R} = -\frac{1}{\lambda} \hat{u}^\mathsf{T}, \qquad \frac{\partial D}{\partial d_R} = -\frac{c}{\lambda} \frac{\partial \dot{\Delta t}_R}{\partial d_R}.\]

For carrier phase:

\[\frac{\partial \Phi}{\partial r_R} = \frac{1}{\lambda}\hat{u}^\mathsf{T}, \qquad \frac{\partial \Phi}{\partial b_R} = \frac{c}{\lambda} \frac{\partial \Delta t_R}{\partial b_R}, \qquad \frac{\partial \Phi_\mathrm{meas}}{\partial N} = 1.\]

Derivatives of transmitter ephemeris interpolation, light-time coupling, Shapiro delay, plasma delay, and C/N0-based noise with respect to receiver state are not included in the current analytic Jacobian.

Online and Precompute Equivalence¶

The online path computes

\[\mathcal{Y}(t_i, x_i) = \operatorname{Compute}(t_i, x_i)\]

one epoch at a time.

The precompute path computes the same measurement epochs for vectors receive_times and user_states:

\[\left\{ \mathcal{Y}(t_i, x_i) \right\}_{i=0}^{M-1}.\]

When a batch ionosphere/plasma provider is configured, precompute first builds all light-time-corrected channels without running the online plasma model:

\[\mathcal{C}_{ij} = \operatorname{BuildChannel}(t_i, x_i, \mathrm{PRN}_j),\]

then calls the batch provider to fill \(\Delta \rho_{\mathrm{plasma},ij}\) for each epoch/channel before evaluating the observables. This keeps visibility, light-time, transmitter clock terms, and batch plasma simulation ordered explicitly, and it avoids doing an online raytrace that would be overwritten by the batch result.

In the staged Lunar GNSS ODTS scenario, this contract is implemented with an explicit file boundary: LunarGnssODTSSimulation::Precompute writes all light-time-corrected links with zero plasma delay, precompute_delays.py fills the delay terms with GCPM/IRI ray tracing, and LunarGnssODTSSimulation::Run merges the delay file back into the truth channels before generating pseudorange, Doppler, and optional TDCP measurements. TDCP is represented as a carrier-range difference in meters and is processed by the UDU stochastic-cloning filter because its measurement model depends on both the current and previous receiver state.