Measurement Simulation

ACONS synthesises observables by combining satellite/user ephemerides, link budgets, and tracking loop noise models. This page documents the signal assumptions and equations implemented in src/measurement_simulator.py, which writes the full measurement catalogue to measurements.json.

Catalogue layout

measurements.json now stores a nested object to keep related quantities together. The top-level keys are:

schema_version: format version (currently 2).
metadata: record count, column list, and scenario hints (e.g., generated measurement types).
measurements: array of per-epoch objects grouping related values.

Each entry contains:

epoch: julian_day / ephemeris_time_seconds.
satellite: satellite id and body-fixed position/velocity vectors.
receiver: body-fixed state plus clock bias/drift.
geometry: visibility flag, azimuth/elevation/zenith, transmit off-boresight, and line-of-sight unit vector.
link_budget: system temperature, $C/N_0$, FSPL, EIRP, and receive gain.
observables: nested range, range_rate, two_way_range, and two_way_range_rate blocks with observed/true values, noise realisations, one-sigma sigmas, two-way delays, and delayed samples (when simulated).
errors: SISE decompositions plus the satellite clock-drift error.

Downstream tools should call measurement_file_manipulation.load_measurement_catalogue(...) to flatten this structure back into the legacy column layout whenever a DataFrame is required.

Mars surface elevation sampling

Rovers or landers that sit on the Martian surface can provide a better initial position by querying the blended HRSC + MOLA DEM stored under data/mars_dem/global/. The helper src.dem.mars_dem.mars_elevation loads only the necessary 2×2 window from Mars_HRSC_MOLA_BlendDEM_Global_200mp_v2.tif and returns either the bilinear-interpolated (default) or nearest-neighbour elevation in metres. Longitudes supplied in [0, 360] are automatically wrapped for rasters that expose [-180, 180] bounds, so callers can stick to whichever convention is convenient.

Download the global HRSC + MOLA blended DEM from https://astrogeology.usgs.gov/search/map/mars_mgs_mola_mex_hrsc_blended_dem_global_200m and place the GeoTIFF where estimation.dem.base_dir/estimation.dem.global_filename point (for example, in configs/scenarios/estimation/ekf_mars.yaml); it is required for src.dem.mars_dem.py.

Unit tests exercise the longitude wrapping and nearest/bilinear sampling using small synthetic GeoTIFF tiles, so changes to the sampling logic remain deterministic.

The default DEM path can be assembled by joining estimation.dem.base_dir and estimation.dem.global_filename when wiring up new scenarios or quick-look scripts.

Reference-point regression checks use the Mars DEM from the estimation config and assert the sampled heights remain within 100 m of the expected values.

For bulk lookups or endpoint-to-endpoint profiles, src/dem/mars_dem.py provides vectorized sampling and simple path generation while retaining nearest/bilinear interpolation.

Terrain projection and slope motion

src/dem/terrain_projection.py adds planet-agnostic BLH/XYZ transforms, terrain projection (project_hold_bl, project_iterative), and slope-constrained motion (step_slope_motion). Coordinate transforms now live in src/coordinates_transformation.py, while the projection/slope utilities remain in src/dem/terrain_projection.py. All BLH inputs use degrees for latitude/longitude. For spherical bodies, pass sphere_radius_m into xyz_to_blh when a mean-radius override is required; otherwise the transform uses the planet parameters supplied to the call. DEM heights are sourced via src.dem.mars_dem.mars_elevation by default, with optional gradient finite differencing.

from src.dem.mars_dem import mars_elevation

height_m = mars_elevation(
    latitude_deg=18.4,
    longitude_deg=77.5,
    dem_path="data/mars_dem/global/Mars_HRSC_MOLA_BlendDEM_Global_200mp_v2.tif",
)

The function returns float("nan") when the coordinate falls outside of the DEM coverage or if the raster reports nodata at the requested location, enabling downstream filters to drop unusable points without flag juggling.

Geometry

For a satellite with position $\mathbf{s}$ and velocity $\mathbf{v}_s$, and a user with position $\mathbf{u}$ and velocity $\mathbf{v}_u$, the simulator evaluates:

\[ \rho = \|\mathbf{s}-\mathbf{u}\|, \qquad \dot\rho = (\mathbf{v}_s-\mathbf{v}_u)^\top\frac{\mathbf{s}-\mathbf{u}}{\rho}, \]

the one-way line-of-sight range and range rate expressed in metres and metres-per-second. These quantities provide the geometric core for the range and Doppler observables stored in measurements.json (fields true_range_m, range_rate_true_mps). The redundant true_range_light_time_m column has been dropped, so downstream tools should rely on true_range_m for the geometric path length.

Signal model

The link budget follows a Friis formulation. For each time sample the simulator evaluates

\[ \frac{C}{N_0} = \frac{\mathrm{EIRP}\cdot L_p \cdot g_r}{k_B\,T_{\mathrm{eq}}}, \]

where $\mathrm{EIRP}$ is the transmit pattern evaluated at the user look angle, $L_p$ is the free-space path loss, $g_r$ is the receive antenna gain, $k_B$ is Boltzmann’s constant, and $T_\mathrm{eq}$ is the equivalent noise temperature computed with Friis’ cascade (antenna plus LNA). measurement.transmitter and measurement.receiver_rf control the gain patterns and RF parameters. The resulting $C/N_0$ (field cn0_dbhz) drives both measurement availability and the thermal noise models described below.

Physical constants (speed of light, Boltzmann constant, ambient temperature) are sourced from configs/environment/constants.yaml via src/constants.py, so update that file to tune the defaults used by the simulator.

Set measurement.receiver_rf.cn0_threshold_dbhz (default 32) to configure the minimum acceptable $C/N_0$. During simulation the tool compares the lowest visible value per satellite against this threshold and emits a warning in simulate.log whenever a spacecraft dips below it, helping highlight geometry/RF issues without editing the code.

Range and Doppler observables

One-way model

Measured pseudoranges can be written as

\[ \hat{\rho}_u = \|\tilde{\mathbf{r}}_s - \mathbf{r}_u\| + c\,(t_u - \tilde{t}_s) + \varepsilon, \qquad \hat{\dot{\rho}}_u = \|\tilde{\dot{\mathbf{r}}}_s - \dot{\mathbf{r}}_u\| + c\,(\dot{t}_u - \tilde{\dot{t}}_s) + \dot{\varepsilon}, \]

where $\tilde{\cdot}$ denotes the transmitted (clock-biased) satellite quantities and $t_u$, $\dot{t}_u$ represent the user clock bias and drift. The simulator realises these observables via

\[ \rho = \|\mathbf{r}_s - \mathbf{r}_u\| + \Delta t_u + \Delta t_{sv} + \varepsilon_{DLL}, $$ $$ \dot{\rho} = \|\dot{\mathbf{r}}_s - \dot{\mathbf{r}}_u\| + \Delta\dot{t}_u + \Delta\dot{t}_{sv} + \varepsilon_{FLL}, \]

with the following components:

$\Delta t_u = b_u(t)$ is the user clock bias in metres. It is obtained by integrating the drift series generated from the Allan deviation curve configured under measurement.oscillator.
$\Delta \dot{t}_u = \dot{b}_u$ is the corresponding clock drift (metres-per-second) converted from the OCXO fractional-frequency noise and applied uniformly across all satellites at a given epoch.
$\Delta t_{sv}$ and $\Delta\dot{t}_{sv}$ model satellite clock/orbit errors (“SISE”). The scenario provides ODTS one-sigma scalars (metres and metres-per-second) under measurement.sise. These sigmas are already defined along the line of sight, so the simulator applies a zero-mean Gaussian draw with that standard deviation to every satellite at a given epoch. Clock bias/drift SISE use the corresponding scalar values provided in the same block and are added on top of the orbit term for one-way observables. The simulator stores the realised series as sise_error_m (combined term for range-like observables), sise_orbit_error_m, sise_clock_bias_error_m, sise_orbit_rate_error_mps, and sat_clock_drift_error_mps. The per-epoch variances are exported alongside the measurements as sise_range_variance_m2, sise_two_way_range_variance_m2, sise_range_rate_variance_mps2, and sise_two_way_range_rate_variance_mps2, enabling the EKF to enlarge its innovation covariance with the same SISE model used by the simulator.
$\varepsilon_{DLL}$ and $\varepsilon_{FLL}$ are the thermal tracking errors derived from the DLL/FLL models detailed below.

Clock and SISE inputs

The simulator expects a single measurement.oscillator entry with Allan deviation samples at two or more averaging times. Values can be supplied at arbitrary $\tau$ and the ingestion stage interpolates across them, ensuring the truth clock matches the specified profile. The measurement.sise block further refines the satellite-side position/velocity and clock errors using the ODTS sigmas, for example:

measurement:
  sise:
    position_sigma_m: 2.89
    velocity_sigma_mps: 7.5e-4
    clock_sigma_m: 7.5
    clock_drift_sigma_mps: 7.5e-4

These values drive both the simulated observables and the innovation covariance used by the EKF.

Two-way observables

Two-way measurements remove the clock terms while doubling the geometric path length. The simulator therefore uses

\[ \rho^{2w} = 2\|\mathbf{r}_s - \mathbf{r}_u\| + 2\Delta t_{sv,p} + c\,t_{cb} + \varepsilon_{TWM}, \qquad \dot{\rho}^{2w} = 2\|\dot{\mathbf{r}}_s - \dot{\mathbf{r}}_u\| + 2\Delta\dot{t}_{sv,p} + \varepsilon_{TWM}, \]

where $t_{cb}$ is the configurable two-way calibration bias (0.5 ns → 0.15 m by default) and $\varepsilon_{TWM} = \sqrt{\varepsilon_{Forward}^2 + \varepsilon_{Return}^2}$, so the two-way noise standard deviation is $\sqrt{2}$ times the one-way value. The ODTS orbit and velocity draws $\Delta t_{sv,p}$, $\Delta\dot{t}_{sv,p}$ are applied on both legs and remain present even after the clock terms cancel. measurements.json records the realised user and satellite clock terms (user_clock_bias_m, user_clock_drift_mps, sise_error_m, sise_orbit_error_m, sise_clock_bias_error_m, sise_range_variance_m2, sise_two_way_range_variance_m2, sise_orbit_rate_error_mps, sise_range_rate_variance_mps2, sise_two_way_range_rate_variance_mps2, sat_clock_drift_error_mps, etc.) alongside the observable values, allowing downstream tools to inspect the individual contributors.

Deep-space assets typically limit uplink time, so the simulator constrains the two-way observables to short windows. measurement.two_way_availability_minutes specifies how long each contact lasts, while measurement.two_way_availability_cadence_minutes defines how often the windows repeat (both default to 60, yielding continuous availability). When the cadence exceeds the duration the simulator gaps out the two-way range and range-rate columns between contacts, ensuring the EKF only ingests the available passes. Setting the duration to 0 removes all two-way data entirely.

To produce a delayed arrival view of two-way data, set measurement.two_way_delay_simulate: true along with measurement.two_way_delay_seconds. The simulator keeps the original two-way columns and adds _delayed companions (for example two_way_range_m_delayed) by interpolating each satellite's time series at t - two_way_delay_seconds (the calibration bias is applied after the delay). The delay value may be fractional. A two_way_delay_seconds column is written to the measurement output and can be edited to inject per-measurement delays. Delayed values are only interpolated within valid measurement segments; outside a valid window the delayed value remains NaN. If you edit two_way_delay_seconds after simulation, the estimate stage rebuilds the _delayed columns when measurement.two_way_delay_simulate: true so they stay consistent.

To run the standard EKF with delayed two-way data (delay not estimated), keep estimation.delayed_twm_enabled: false. When delayed two-way columns are present, the estimate stage swaps the _delayed two-way columns into the base two-way fields before filtering.

measurement:
  types: [range, range_rate, two_way_range]
  two_way_availability_minutes: 10     # 10-minute windows
  two_way_availability_cadence_minutes: 30  # repeating every 30 minutes
  two_way_selection_strategy: window_locked  # lock to the highest-C/N0 link at window start

measurement.two_way_selection_strategy controls whether the highest-$C/N_0$ satellite is chosen at every epoch inside a contact (per_epoch, legacy behaviour) or only once when the window opens (window_locked, keeping the same spacecraft throughout the window). The latter better emulates ground stations that avoid reconfiguring beams multiple times within the same pass.

Thermal noise and tracking loops

The measurement block in the scenario file exposes the delay-lock (DLL) and frequency-lock (FLL) loop parameters. The simulator uses the Methodology jitter model to convert $C/N_0$ into measurement noise:

\[ \sigma_\text{DLL} = \frac{c/f_\text{code}}{2\,\Delta} \sqrt{\frac{B_L}{C/N_0}\left(1+\frac{1}{T_i\,C/N_0}\right)}, \]

\[ \sigma_\text{FLL} = \frac{\lambda}{2T_i} \sqrt{F\frac{B_L}{C/N_0}+\frac{1}{T_i\,(C/N_0)^2}}, \]

where

$B_L$ — loop bandwidth (0.5 Hz by default),
$T_i$ — coherent integration time (20 ms),
$\Delta$ — early–late spacing (1 chip),
$\lambda = c/f_c$ — carrier wavelength,
$F$ — PLL/FLL scaling factor (1 above 35 dB-Hz, 2 below, configurable).

The resulting $\sigma_\text{DLL}$ (metres) and $\sigma_\text{FLL}$ (converted to m/s) populate the range_noise_std_m and range_rate_noise_std_mps fields. The EKF reads these values directly from measurements.json.

Measurement diagnostics

measurement_simulator.py retains the full set of metadata per observation (satellite/user position and velocity, LOS vector, azimuth/elevation, visibility flag). The reporting layer applies a robust Median Absolute Deviation (MAD) filter before generating summaries and plots, ensuring extreme outliers do not skew the diagnostics. See Outputs & Reporting for the derived CSV/plot artefacts.

Trace logging

Set --log-level TRACE when running acons simulate … to include exemplar measurement breakdowns in simulate/simulate.log. The simulator logs one representative one-way and two-way range sample plus the corresponding range-rate observables, highlighting the geometric term, user clock bias or drift, satellite orbit/clock errors, calibration bias, and DLL/FLL thermal noise that sum to each observation.

Implementation notes

The simulator code mirrors the structure described above. Geometry extraction, link-budget evaluation, and observable synthesis live in dedicated helpers (_satellite_geometry, _link_budget_metrics, _generate_range_observables, _generate_rate_observables), while simulate_measurements orchestrates the per-satellite loop. The split keeps the numeric models focused and makes it easier to extend the simulator with alternative antennas, additional measurement types, or unit tests that exercise a single stage of the pipeline.