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.

Trace logging

Run poetry run acons simulate --log-level TRACE … to enrich simulate/simulate.log with structured entries that expand on the INFO-level events. TRACE now includes the resolved receiver longitude/latitude/height/source and a bullet-style summary of the measurement configuration (observable types, carrier/chip rates, baseline code noise, two-way availability/delay/selection strategy, RNG seed). These fields mirror the scenario YAML so you can verify the simulator is using the intended parameters before generating measurements.json.

Planetary 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.planet_dem.planet_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). planet_elevation will build the geodetic CRS from the scenario’s celestial_body (or from explicit body axes), so projected DEMs stay aligned with the configured target body. If the GeoTIFF CRS already includes body axes (+a/+b, +R, or a WKT SPHEROID), those values are preferred for full CRS fidelity.

For quick one-off checks, use the standalone sampler:

poetry run python scripts/dem_sample.py data/mars_dem/global/Mars_HRSC_MOLA_BlendDEM_Global_200mp_v2.tif 24.6 50.4

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/planet_dem.py provides vectorized sampling and simple path generation while retaining nearest/bilinear interpolation.

Altimeter measurements (EDL only)

When measurement.altimeter is enabled for an EDL scenario, the simulator evaluates the altimeter model in src/sensors/altimeter.py and records the results in the measurement catalogue. The output columns include altimeter_agl_m (measured altitude above ground level), altimeter_true_agl_m (noise-free AGL), altimeter_noise_std_m (per-epoch sensor sigma used by the EKF), altimeter_terrain_h_m (DEM terrain height used by the sensor), altimeter_vehicle_height_datum_m (vehicle height above the DEM datum after undulation correction), altimeter_undulation_m (body-dependent areoid/selenoid undulation), and altimeter_valid (boolean validity flag). The simulator also persists individual sensor state flags as altimeter_flag_out_of_range, altimeter_flag_dropout, altimeter_flag_hold, and altimeter_flag_latency_wait. Measurements are gated by measurement.altimeter.enable_below_m (default 1000 m), so values above the activation height are left as NaN. For EDL users the simulator evaluates the sensor on the time-varying receiver trajectory reconstructed from receiver_x_m/y_m/z_m. The body-fixed XYZ state is converted to planetocentric latitude/longitude and a body-dependent reference-surface height, then corrected by a planetary undulation model before DEM sampling. When a gravity-model provider is available, the reference-surface height is formed with that model's native reference radius (r0); the scenario radius remains a fallback:

h_vehicle_datum = h_reference - N

AGL = h_vehicle_datum - h_terrain_datum

where N is the planetary undulation obtained from src/dem/planetary_undulation.py. The terrain height used to form AGL comes from the DEM referenced by the active estimation config. Altimeter noise/bias states are seeded with the same measurement.rng_seed used by the rest of the simulation pipeline so repeated runs stay deterministic.

Undulation providers are selected automatically from celestial_body:

Mars uses pyshtools.datasets.Mars.GMM3()
Moon uses pyshtools.datasets.Moon.GRGM1200B()

The pyshtools datasets are loaded lazily and cached in the user cache directory after first use.

Enable the altimeter in a measurement YAML (EDL only):

measurement:
  altimeter: !include components/altimeter.yaml

EDL scenario example:

configs/scenarios/mars_edl.yaml combines the EDL user definition, altimeter-enabled measurement model, and EDL DEM settings.

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.planet_dem.planet_elevation by default, with optional gradient finite differencing.

from src.dem.planet_dem import planet_elevation

height_m = planet_elevation(
    latitude_deg=18.4,
    longitude_deg=77.5,
    dem_path="data/mars_dem/global/Mars_HRSC_MOLA_BlendDEM_Global_200mp_v2.tif",
    scenario_path="configs/scenarios/estimation/ekf_mars.yaml",
)

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.

Rover trajectory-driven receiver state

For user.type: rover with user.vehicle_trajectory_path set, the simulator now samples the receiver state at every epoch from the loaded vehicle ephemeris (including GeoJSON trajectories like data/vehicle/RoverDrivePath.geojson). The trajectory is interpolated on the simulation ET axis and written into receiver_x_m, receiver_y_m, receiver_z_m, receiver_vx_mps, receiver_vy_mps, and receiver_vz_mps per epoch.

Vehicle trajectories are transformed into the scenario frame defined by frames.body_fixed before interpolation. This means STK *.e files declared in ICRF/J2000 are automatically rotated into the active body-fixed frame when the required SPICE kernels are furnished.

This means range/range-rate observables are generated against the moving rover state rather than a single fixed receiver location. For user.type: rover and user.type: edl, the initial lon/lat/height are also derived directly from user.vehicle_trajectory_path; user.location is not used in these moving-user modes.

The scenario time window is authoritative for moving users. Absolute timestamps embedded in user.vehicle_trajectory_path are ignored during simulation; the first trajectory sample is aligned to the scenario start, and only the elapsed timing between support points is preserved. This allows trajectory files such as edl_insight.csv to be reused with a different scenario epoch without changing the physical shape of the motion profile.

To inspect the interpolation quality itself, run py -3 scripts/analyze_vehicle_interpolation.py. The diagnostic compares a sparse rover path against a denser reference trajectory, plots truth vs. interpolated x/y/z, and writes the resulting figure and error tables to outputs/vehicle_interpolation/.

When the vehicle trajectory is too sparse relative to the requested simulation cadence, the simulator now inserts additional support points only between existing samples before evaluating the rover state on the observation grid. This enables observation generation at the configured cadence without extrapolating outside the real trajectory span. The inserted support points are logged as artificial and the effective rover velocity used for simulation is now derived from the interpolated position timeline, so receiver_v* remains consistent with receiver_x/y/z.

The simulator stores the trajectory preparation summary in measurements.json.metadata.trajectory_meta and restores it when reloading the catalogue through measurement_file_manipulation.load_measurement_catalogue(...).

For user.type: rover, the support trajectory is additionally projected onto the DEM surface before interpolation onto simulation epochs. This enforces reference-surface consistency between rover truth generation and DEM-based estimation constraints. The projection status and height-shift statistics are recorded in trajectory_meta under dem_projection_* keys. DEM projection sampling is executed in vectorized mode (batch DEM lookup over all support points) to reduce simulation runtime; trajectory_meta.dem_projection_sampling_mode reports the active mode.

For Mars rover studies, configs/scenarios/mars_rover_case.yaml provides a dedicated scenario entry point that keeps the rover-specific user include explicit. The simulator logs the warning below whenever such artificial support points are created:

WARNING: Trajectory sampling is too sparse. Additional points were generated via interpolation.
Generated trajectory should be treated as artificial.

Trajectory handling now follows the overlap between the requested simulation window and the vehicle ephemeris:

if the trajectory is longer than the simulation window, it is trimmed to the requested range
if the trajectory is shorter than the simulation window, the effective simulation ends at the trajectory end instead of holding the last position constant
if the overlap leaves fewer than 5 effective simulation epochs on the scenario time grid, the simulator now stops with a configuration error instead of running a near-empty moving-user case

Simulation plotting adds simulate/plots/receiver_trajectory_sampling.png, which distinguishes the original support points from any interpolated support points and marks the trajectory end inside the effective simulation window.

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. The simulator builds a fractional-frequency series $y(t)$ by combining white FM, flicker FM, and random-walk FM basis noises, then fits their amplitudes (via non-negative least squares) so the Allan deviation matches the measurement.oscillator targets. The time bias $x(t)$ is obtained by integrating $y(t)$ and then scaled by $c$ to express the bias in metres.
$\Delta \dot{t}_u = \dot{b}_u$ is the corresponding clock drift (metres-per-second) obtained by scaling the simulated fractional-frequency series by $c$ 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. For each satellite the simulator draws a normally distributed constant bias $b$ for the position, velocity, clock bias, and clock drift terms. These satellite-dependent biases use fixed sigmas of 0.01 m (position), 0.1 mm/s (velocity), 0.01 m (clock bias), and 0.1 mm/s (clock drift). For example $\Delta r_{sv} = \mathcal{N}\big(b_{sv,r}, \sigma_{sv,r}\big)$ with $b_{sv,r} = \mathcal{N}(0, 0.01\,\text{m})$, and the per-epoch samples share the same bias. The vector position/velocity errors are added to every component of the satellite state before projecting onto the line of sight, ensuring the range and Doppler observables pick up the correct contribution from each axis. 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 simulator fits the white/flicker/random-walk FM amplitudes so the realised $y(t)$ matches those targets over the available time span. Allan deviation fitting uses allantools, while flicker FM synthesis relies on colorednoise. 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 reuse the same geometric path length as the one-way model but cancel the user clock terms. The simulator therefore uses

\[ \rho^{2w} = \|\mathbf{r}_s - \mathbf{r}_u\| + \Delta t_{sv,p} + c\,t_{cb} + \varepsilon_{TWM}, \qquad \dot{\rho}^{2w} = \|\dot{\mathbf{r}}_s - \dot{\mathbf{r}}_u\| + \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 user clock cancels, while the satellite clock drift term is excluded for two-way range-rate. 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.

Error-budget plots

Each simulation run also writes simulate/plots/error_budget.png beneath the configured --run-dir. The figure contains one subplot per enabled observable (range, range_rate, two_way_range, two_way_range_rate) and stacks the mean absolute contribution of each modeled error source per satellite: satellite orbit and clock SISE draws, the DLL/FLL thermal noise, and the injected two-way calibration bias. User clock bias and drift are excluded from this plot. This provides a quick visual budget of which satellites or error sources dominate each observable without digging through measurements.json.

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 per-observable noise statistics in simulate/simulate.log. Instead of a single sample, the simulator now reports (for each enabled measurement type) the sample count, thermal-noise RMS/max, SISE orbit/clock RMS/max, the distribution of DLL/FLL-derived σ values (mean/min/max), and the nominal SISE σ inputs pulled from the scenario. Two-way ranges also note the injected calibration bias. These summaries quantify how the configured error sources propagate into measurements.json, making it easy to validate that different observables are picking up the expected error budgets. Each satellite that clears the C/N₀ threshold also receives its own TRACE block (emitted next to the INFO Generated measurements line) so you can audit per-spacecraft observation mixes without waiting for the end-of-run aggregate summary.

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.