Forward-Backward Map Matching for Airborne Magnetic Anomaly Navigation

Joseph O'Leary; Tapio Simula; Tom Kimpson; Corey Davey; Wenchao Li; Andrew D. Greentree; Allison Kealy

eess.SPOct 27, 2026

Forward-Backward Map Matching for Airborne Magnetic Anomaly Navigation

Authors: Joseph O'Leary, Tapio Simula, Tom Kimpson, Corey Davey, Wenchao Li, Andrew D. Greentree, Allison Kealy

Organizations: Innovative Planet Research Institute, Swinburne University of Technology, Melbourne, Australia · School of Physics, University of Melbourne, Melbourne, Australia · Optical Sciences Centre, Swinburne University of Technology, Melbourne, Australia · School of Mathematics and Statistics, The University of Melbourne, Melbourne, Australia ·

Abstract

Airborne magnetic anomaly navigation provides a passive, unjammable alternative to Global Navigation Satellite Systems positioning by matching magnetometer measurements against a lithospheric magnetic anomaly map. Existing map-matching algorithms either assume a Gaussian posterior (extended Kalman filter), require Monte Carlo sampling (particle filter), or discard posterior uncertainty in favour of a single optimal path (Viterbi algorithm). Here we formulate map matching as inference on a chain-structured factor graph and solve it segment by segment using the forward-backward algorithm, which yields a discrete posterior over candidate positions and its Shannon entropy $H_t$ at time $t$ . Over the Renfrew\_395 map, $10^3$ synthetic trials yield a median position error of $32.7$ m, $91\%$ below the median inertial drift of $362$ m. On these trials and $10^3$ flight paths, the median error is $43\%$ below the extended Kalman filter and particle filter. Across these paths, a threshold on segment-level entropy detects $47\%$ of segments with errors above two map cells and $94\%$ of those above three, at a $10\%$ false-alarm rate. The foregoing results indicate that $H_t$ provides a truth-free diagnostic of map-conditioned navigation ambiguity.

Figures & tables

Figure 1: A chain-structured factor graph with prior

f_{0}

, measurement

f_{m}^{(\ell)}

, and transition

f_{\tau}

factors. Hidden states

s_{t_{1}},\ldots,s_{t_{T}}

(solid circles) are connected by transition factors

f_{\tau}

(filled squares). Two families of measurement factors,

f_{m}^{(1)}

and

f_{m}^{(2)}

(dashed open squares), link each hidden state to corresponding observations

z_{t}^{(1)}

and

z_{t}^{(2)}

(dashed circles). A prior factor

f_{0}

(open square with solid border) is attached to the initial state

Figure 2: Forward and backward sweeps on the chain of Figure 1 , with one measurement factor per timestep. The forward variable

\alpha_{t_{k}}

collects the prior factor

f_{0}

and the factors up to

t_{k}

(blue), and the backward variable

\beta_{t_{k}}

collects the factors after

t_{k}

(orange). The forward sweep computes

\alpha_{t_{1}},\dots,\alpha_{t_{k}}

by Equations ( 4 ) and ( 5 ), and the backward sweep computes

\beta_{t_{T}},\dots,\beta_{t_{k}}

by Equations ( 6 ) and ( 7 ). The normalised product of

\alpha_{t_{k}}

and

\beta_{t_{k}}

s_{t_{k}}

(shaded) is the marginal posterior

\gamma_{t_{k}}

of Equation ( 3 )

Figure 3: Schematic of the candidate set $\mathcal{S}_{t_{k}}$ on a $5\times 5$ grid at three successive timesteps $t_{1}$ , $t_{2}$ , and $t_{3}$ . Each grid is centred on the map node nearest to $\hat{\bm{s}}^{\rm INS}_{t_{k}}+\Delta\hat{\bm{s}}$ (grey points, dashed curve), the INS position plus the correction carried from the previous segment [Equations ( 12 ) and ( 13 )]. The shading is the marginal posterior density $p(\bm{s}_{t_{k}}\mid z_{\mathcal{T}})$ , and $\gamma_{t_{k}}(j)$ [Equation ( 3 )] is its value at the cell centres, normalised over the grid. The marginal MAP estimate $\hat{\bm{s}}^{\rm MAP}_{t_{k}}$ [Equation ( 9 ); blue squares, solid curve] is the cell of highest $\gamma_{t_{k}}$ at each timestep. The Viterbi path of Appendix A , the most probable sequence of cells over the whole segment, does not in general pass through these cells. The MMSE estimate $\hat{\bm{s}}^{\rm{MMSE}}_{t_{k}}$ [Equation ( 8 ); orange points, solid curve] is the posterior-weighted average over the candidate set. The $5\times 5$ grid is presented as a visual aid; the results in Sections 4 and 5 adopt $n=101$

Figure 4: The Renfrew_395 aeromagnetic anomaly map $B_{\rm anom}(\phi,\lambda)$ (units: nT), overplotted with the ground-truth trajectory $\bm{s}_{t_{k}}$ (black, solid curve) and the INS position estimate $\hat{\bm{s}}_{t_{k}}^{\rm INS}$ (red, dashed curve). The start and finish positions are marked with a white circle and yellow diamond, respectively. Geodetic longitude $\lambda$ (units: degrees) is reported on the horizontal axis and geodetic latitude $\phi$ (units: degrees) on the vertical axis. The colour bar spans the anomaly field at the flight altitude $h=1000$ m over the plotted region, from its smallest to its largest map value

Parameter	Symbol	Value	Units
Initial position uncertainty	$\sigma_{\rm pos,0}$	100	m
Initial velocity uncertainty	$\sigma_{\rm vel,0}$	0.01	m s -1
Initial attitude uncertainty	$\sigma_{\rm att,0}$	0.01	deg
Velocity random walk	$\sigma_{\rm VRW}$	$2.38\times 10^{-4}$	m s -3/2
Angular random walk	$\sigma_{\rm ARW}$	$5.81\times 10^{-7}$	rad s -1/2
Accelerometer bias	$\sigma_{\rm a}$	$2.45\times 10^{-4}$	m s -2

Table 1: INS error model parameters for the representative realisation.

Parameter	Symbol	Value	Units
White measurement noise variance	$\sigma_{\rm meas}^{2}$	1.0	nT 2
FOGM amplitude	$\sigma_{S}$	1.0	nT
FOGM correlation time	$\tau_{S}$	600	s

Table 2: Magnetic measurement model parameters.

Parameter	Symbol	Value	Units
Measurement noise (evidence maximum)	$\sigma_{\rm mag}$	$0.82$	nT
Velocity uncertainty (evidence maximum)	$\sigma_{v}$	$3.26$	m s -1
Grid dimension	$n$	$101$	–
Map grid spacing	$\Delta_{\rm cell}$	$29.3$	m
Segment length (derived)	$T$	$40$	–
Sampling interval	$\Delta t_{\rm FG}$	5.0	s

Table 3: Factor graph parameters for the representative test case.

Figure 5: Position error

e_{\rm pos}(t)

(units: m) versus time

t

(units: min) for the representative test case of Section 4.2 . The error of the INS position estimate is plotted as a red, dashed curve and that of the forward-backward MMSE estimate (Equation 8 ) as a cyan, solid curve. The errors of the forward-only PMF (orange, solid curve), the bootstrap particle filter with

10^{5}

particles (green, dotted curve), and the Viterbi path (purple, solid curve) of Appendix A are overplotted for comparison. The

M=N_{\rm FG}/T=6

segment boundaries (Section 3.7 ) are overplotted as grey, dashed, vertical lines. The vertical axis is plotted on a logarithmic scale and the horizontal axis on a linear scale. The factor graph parameters adopted in this section are listed in Table 3

Figure 6: Shannon posterior entropy

H_{t}

[Equation ( 10 ); units: nats] for the representative test case of Section 4.2 . Left panel:

H_{t}

versus the per-timestep Fisher trace

\operatorname{tr}(J_{t})

(units: km -2 ) on a logarithmic horizontal axis, with one point per timestep (

N_{\rm FG}=240

). Middle panel: segment-level median entropy

H^{(m)}_{\rm med}=\operatorname{median}_{t\in\mathcal{T}^{(m)}}H_{t}

versus the segment-level minimum eigenvalue of the accumulated Fisher information

\lambda_{\min}[J^{(m)}]

(units: km -2 ) on a logarithmic horizontal axis, with one point per segment (

M=6

). The vertical bars denote the 16th to 84th percentile of

H_{t}

within each segment. Right panel:

H_{t}

versus position error

e_{\rm pos}(t)

(units: m) on a linear horizontal axis, with one point per timestep (

N_{\rm FG}=240

)

Figure 7: Per-trial median posterior entropy

H_{\rm med}=\operatorname{median}_{t}H_{t}

[Equation ( 10 ); units: nats] versus per-trial median position error

e_{\rm med}=\operatorname{median}_{t}e_{{\rm pos},t}

(units: m) for

10^{3}

independent realisations of the INS and magnetometer error processes on a fixed flight path (Section 4.2 ). Each point corresponds to a single sensor-error realisation, with the medians taken over all

N_{\rm FG}=240

timesteps within each trial. Points are coloured by percentile bands of

e_{\rm med}

across trials: green points denote trials below the 84th percentile (

e_{\rm med}\leq 40.9

m), yellow points denote trials between the 84th and 97.5th percentiles (

40.9

–

50.6

m), orange points denote trials between the 97.5th and 99.85th percentiles (

50.6

–

64.6

m), and red points denote trials above the 99.85th percentile (

>64.6

m). Both axes are plotted on linear scales. The factor graph configuration (

n=101

T=40

M=6

) is identical across all trials

Figure 8: Within-segment median position error $e^{(m)}_{\rm med}=\operatorname{median}_{t\in\mathcal{T}^{(m)}}e_{{\rm pos},t}$ (units: m) for $10^{3}$ distinct flight paths over the Renfrew_395 map with a common sensor-error realisation. Each trajectory contributes $M=6$ segments, yielding $6\times 10^{3}$ segment-level samples in total. In both panels, individual segment observations are overlaid as faint markers, the solid black curve traces the binwise conditional median, and the shaded band spans the 16th to 84th percentile range within each bin. Left panel: the horizontal axis is the segment-level minimum eigenvalue of the accumulated Fisher information $\lambda_{\min}[J^{(m)}]$ (units: km -2 ), plotted on a logarithmic scale. Right panel: the horizontal axis is the within-segment median posterior entropy $H^{(m)}_{\rm med}=\operatorname{median}_{t\in\mathcal{T}^{(m)}}H_{t}$ [Equation ( 10 ); units: nats], plotted on a linear scale. The factor graph configuration ( $n=101$ , $T=40$ , $M=6$ ) is identical across all trials; see Table 3

	Median error (m)
Estimator	Sensor-error realisations	Flight paths	Run time (s)
Factor graph, MMSE	32.8	23.8	0.022–0.032
Factor graph, marginal MAP	35.5	27.0	0.022–0.032
Viterbi path	38.3	30.5	0.019–0.022
PMF (forward only)	57.7	41.9	0.013–0.020
EKF	57.9	41.6	$<10^{-4}$

Table 4: Median position error of each estimator, pooled over all timesteps and trials, for the

10^{3}

sensor-error realisations of Section 5.1 and the

10^{3}

flight paths of Section 5.2 , and run time per timestep on one core. The EKF and the particle filter use the same two-dimensional model, parameters, and measurements as the factor graph. The PMF, the Viterbi path, and the factor graph use the same grids and pruning.

Figure 9: Detection rate

P_{\rm D}(\kappa)

versus false-alarm rate

P_{\rm FA}(\kappa)

of a threshold

\kappa

on the within-segment median posterior entropy

H^{(m)}_{\rm med}=\operatorname{median}_{t\in\mathcal{T}^{(m)}}H_{t}

[Equation ( 10 ); units: nats], for the

6\times 10^{3}

segments of the

10^{3}

sensor-error realisations of Section 5.1 (left panel) and the

10^{3}

flight paths of Section 5.2 (right panel). A segment has a large error when its median position error

e^{(m)}_{\rm med}

exceeds

40.9

m (orange curves), two cells, i.e.

58.5

m (blue curves), or three cells, i.e.

87.8

m (black curves). Solid curves: entropy of the posterior

\gamma_{t_{k}}

of the factor graph. Dotted curves: entropy of the normalised forward variable

\alpha_{t_{k}}

. The grey, dashed diagonal corresponds to a flag independent of the error

Appendix figures & tables3 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 10: Effect of the pruning threshold

\alpha_{\rm prune}

10^{2}

sensor-error realisations of the representative test case of Section 4 . Left panel: exclusion rate

P_{\rm excl}

, the percentage of timesteps at which the cell nearest the true position is excluded from the active set

\mathcal{A}_{t_{k}}

. Middle panel: entropy bias

\Delta H

(units: nats), the change in the flight-mean entropy relative to the unpruned run. Right panel: change in the median position error

\Delta e_{\rm pos}

(units: m) relative to the unpruned run. Points mark the median over the realisations, and bars span the 16th to 84th percentiles. The dashed vertical line marks

\alpha_{\rm prune}=10^{-2}

, the setting adopted in this paper. The horizontal axes are plotted on a logarithmic scale

Figure 11: Sensitivity of the median position error $e_{\rm pos}$ (units: m) to the noise parameters for the representative test case of Section 4 . Left panel: median error versus $\sigma_{v}$ (units: m s -1 ) at $\hat{\sigma}_{\rm mag}=0.82$ nT. Middle panel: median error versus $\sigma_{\rm mag}$ (units: nT) at $\hat{\sigma}_{v}=3.26$ m s -1 . Right panel: logarithm of the evidence of Equation ( 37 ) (units: nats) on the $5\times 5$ grid of Appendix C.2 , relative to its maximum; values below $-100$ nats are shown at $-100$ nats. The orange circle marks the evidence maximum in each panel. The left and middle panels share the vertical axis. The axes carrying $\sigma_{v}$ and $\sigma_{\rm mag}$ are plotted on a logarithmic scale

Case	$\bm{\delta B_{\rm rms}}$ (nT)	Median error (m)	Change (m)
No mismatch	0	33.2	–
Measurements from Renfrew_555	0.11	33.4	$0.0$ ( $-1.5$ , $+1.7$ )
Map continued to $h+50$ m	1.96	80.9	$+45.7$ ( $+31.8$ , $+69.9$ )
Map continued to $h+100$ m	3.82	128.5	$+90.9$ ( $+73.7$ , $+132.4$ )

Appendix

Table 5: Magnetic map mismatch on the first

10^{2}

trials of Section 5.1 , each paired with the run without mismatch.

\delta B_{\rm rms}

is the root-mean-square map error along the true track at the flight altitude. The median error is pooled over all timesteps and trials. The change is the per-trial median error minus that of the paired run: median over trials, with the 16th and 84th percentiles in brackets. The median INS error on these trials is 341 m.

Map-Based Localization

Map-Based Localization estimates a device or robot’s position and orientation by matching sensor observations to stored maps.

See all 107 papers

Airborne Magnetic Anomaly Navigation with Neural-Network-Augmented Online Calibration
Mar 9, 2026Antonia Hager, Sven Nebendahl, Alexej Klushyn +3Online Calibration GNSS-Denied Localization
Airbus, 82024 Taufkirchen, Germany. · Norwegian University of Science and Technology, 7034 Trondheim, Norway.
Bio-Inspired Magnetic Navigation Using Gradient Descent With Rao-Blackwellized Particle Filtering
Oct 16, 2026Zayed Herma, Andrea Vassilev, Philippe Cardin +1Robot Navigation GNSS-Denied Localization
Univ. Grenoble Alpes, CEA, Leti, 38000 Grenoble, France · Univ. Grenoble Alpes, CNRS, ISTerre, 38000 Grenoble, France
Observation Quality in DSM-Assisted Aerial Localization: A Benchmark and Base-Anchored Uncertainty Calibration
Oct 13, 2026Pengyu Yin, Jie Pan, Junxiang Tan +3Uncertainty Calibration Map-Based Localization
Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China · School of Electronic Engineering, Beijing University of Posts and Telecommunications, Beijing 100876, China

Mar 9, 2026J/K move · Enter open · S save

Airborne Magnetic Anomaly Navigation with Neural-Network-Augmented Online Calibration

Antonia Hager, Sven Nebendahl, Alexej Klushyn, Jasper Krauser+2

Airbus, 82024 Taufkirchen, Germany. · Norwegian University of Science and Technology, 7034 Trondheim, Norway.

Online Calibration GNSS-Denied Localization Model Calibration Extended Kalman Filtering

Open paper AskPDF arXiv

Airborne Magnetic Anomaly Navigation (MagNav) provides a jamming-resistant and robust alternative to satellite navigation but requires the real-time compensation of the aircraft platform's large and dynamic magnetic interference. State-of-the-art solutions often rely on extensive offline calibration flights or pre-training, creating a logistical barrier to operational deployment. We present a fully adaptive MagNav architecture featuring a "cold-start" capability that identifies and compensates for the aircraft's magnetic signature entirely in-flight. The proposed method utilizes an extended Kalman filter with an augmented state vector that simultaneously estimates the aircraft's kinematic states as well as the coefficients of the physics-based Tolles-Lawson calibration model and the parameters of a Neural Network to model aircraft interferences. The Kalman filter update is mathematically equivalent to an online Natural Gradient descent, integrating superior convergence and data efficiency of state-of-the-art second-order optimization directly into the navigation filter. To enhance operational robustness, the neural network is constrained to a residual learning role, modeling only the nonlinearities uncorrected by the explainable physics-based calibration baseline. Validated on the MagNav Challenge dataset, our framework effectively bounds inertial drift using a magnetometer-only feature set. The results demonstrate navigation accuracy comparable to state-of-the-art models trained offline, without requiring prior calibration flights or dedicated maneuvers.

Abstract

Map-Based Localization

Explore similar work

Airborne Magnetic Anomaly Navigation with Neural-Network-Augmented Online Calibration

Bio-Inspired Magnetic Navigation Using Gradient Descent With Rao-Blackwellized Particle Filtering

Observation Quality in DSM-Assisted Aerial Localization: A Benchmark and Base-Anchored Uncertainty Calibration