A Horizontal Accuracy Metric for Magnetic Navigation

Prasenjit Sengupta

doi:10.33012/navi.717

PDF

Abstract

This paper develops a navigation accuracy metric for magnetic navigation when adapted to civil aviation. Metrics that are currently used may not directly apply to magnetic navigation because the assumptions behind these metrics are based on the statistics of more conventional navigation modalities where the uncertainty distribution in the horizontal plane is typically radially symmetric. Magnetic navigation challenges these assumptions. New standard deviation limits based on the probability of exceedance of radial error are derived. Use is made of quadratic forms over random variables, specifically the Hoyt distribution and its associated probability density function and cumulative density function. Although the density functions lack closed-form solutions, new approximations for this distribution are utilized in order to make rapid computation possible, enabling their use for navigation purposes. When used in conjunction with robust numerical methods, the new approach accurately calculates the bounds on error distributions that would meet the requirements associated with Performance Based Navigation.

Keywords

1 INTRODUCTION

Magnetic Navigation (MagNav) is a rising form of alternative navigation which uses a map of Earth’s crustal magnetic field over a given area of interest and compares the mapped value at an initial estimate of the position to the measured output of a quantum magnetometer in order to update the estimate.

Earth’s crustal magnetic field (or anomaly field) results from permanent or induced magnetization in the crust. It has a magnitude of approximately 100 nT, with spatial variation of up to 200 nT at an altitude of approximately 1 km (Gnadt, 2022a). While the crustal field is much smaller in magnitude than the Earth’s core field (approximately 50,000 nT), its spatial variability and relative temporal stability (approximately 1 nT per year) make it suitable for navigation using sensitive magnetic scalar and vector sensors. Canciani and Raquet (2016) provide an extensive background and history of the MagNav problem and show that a navigation accuracy of up to 13 m can be achieved using a marginalized particle filter, high-quality magnetic anomaly map, and low flight altitudes. Historically, MagNav has been developed for defense applications (Canciani, 2022), but there is growing interest in using it for civil aviation to protect against jamming and spoofing of global navigation satellite systems (GNSS).

For civil aviation applications, navigation accuracy requirements are outlined by the Area Navigation (RNAV) and Required Navigation Performance (RNP) specifications, which are part of the Performance-Based Navigation (PBN) framework (Federal Aviation Administration, 2024). Specifically, the Federal Aviation Administration (2024) definitions state that, “for both RNP and RNAV designations, the numerical designation refers to the lateral navigation accuracy in nautical miles, which is expected to be achieved at least 95 percent of the flight time by the population of aircraft operating within the airspace, route, or procedure.” In the present work, RNAV and RNP are used interchangeably to describe navigation accuracy, although RNP also includes requirements for integrity and continuity, which are beyond the scope of this discussion.

Functionally, meeting this standard in civil aviation settings translates to showing that the uncertainty in the lateral or cross-track error, expressed as a standard deviation, remains below a given threshold. In the civil aviation context, the navigation system uses a filter to estimate an aircraft’s position, using either signals-in-space (e.g., the Global Positioning System (GPS) (SC-159, 2020a, 2020b)) or a combination of navigation aids (NAVAIDs) (e.g., very high frequency omnidirectional range and distance measuring equipment) (Federal Aviation Administration, 2024).

According to “ICAO Doc 9613, Performance-based Navigation (PBN) Manual” (2023), there are currently no formal specifications for longitudinal error apart from Navigation System Error (Section 2.2.2) and RNP 4 (Part C). However, Section 1.6 of the manual discusses the potential inclusion of longitudinal accuracy requirements in future specifications, and Section 2.3.3.3 notes that monitoring and alerting systems vary in how they handle these errors: some apply the limit to the standard deviation of each axis (cross-track and along-track) separately, whereas others apply the 95 percent limit to the combined (2D) radial error. This distinction is nuanced, as the error distribution along individual axes can be assumed to be Gaussian, while radial errors are non-negative and therefore not Gaussian. However, this distinction generally does not affect GPS-based or NAVAID-based positioning because the uncertainty introduced by both approaches does not depend on the direction of travel except in high-latitudes regions (for GPS) and areas with poor NAVAID coverage.

GPS-based navigation systems explicitly account for this lack of directional sensitivity. For example, GPS Dilution of Precision (DOP) (Langley, 1999) evaluates vertical accuracy (VDOP) and horizontal accuracy (HDOP) separately, but HDOP is usually not decomposed into the north-east or longitudinal-lateral components unless a specific application requires it. For example, Satellite-Based Augmentation Systems (SBAS) use the major principal axis variance when calculating the horizontal protection level (HPL) for the position solution (SC-159, 2020a). Integrity monitoring using the Advanced Receiver Autonomous Integrity Monitoring (ARAIM) algorithm (Blanch et al., 2015) similarly involves calculating the HPL using the norm of the HPL along the two horizontal axes, an approach that implicitly assumes that both axes have the same error distribution.

In contrast, MagNav systems are sensitive to the effect of direction on positioning accuracy. This becomes evident when considering the spatial variation in magnetic anomaly values across different regions of the Earth. Figure 4 shows magnetic anomaly measurements sourced from the North American Magnetic Anomaly Group (2002) (NAMAM) at altitudes of 1000 m and 4000 m relative to the ellipsoidal Earth. At higher altitudes, the measurable magnetic anomaly generally has smaller magnitudes and gradients, leading to lower MagNav positioning accuracy. However, strong directionality in the magnetic anomaly gradient, regardless of its absolute magnitude, can occur at any altitude.

This effect is further illustrated in Figures 1 through 3, which show the magnetic anomaly gradient along three different heading angles with respect to true north: 0, 60, and 120 degrees. In these figures, the gradient at each point is calculated with respect to the latitude and longitude at that point and then projected onto the unit vector along the desired heading. Regions with stronger gradients typically have better MagNav performance (i.e., accuracy). In regions with a low gradient magnitude, the magnetic anomaly does not change significantly between two measurement updates, so the estimator cannot infer how the current estimate should be updated. This classic observability problem is encountered in any filtering application (Crassidis & Junkins, 2004).

FIGURE 1

Magnetic anomaly and Cramer-Rao lower bound, heading = 0 deg

FIGURE 2

Magnetic anomaly and Cramer-Rao lower bound, heading = 60 deg

FIGURE 3

Magnetic anomaly and Cramer-Rao lower bound, heading = 120 deg

FIGURE 4

Magnetic anomaly at 1000 m (left) and 4000 m (right)

Figures 1 through 3 additionally show how the magnetic anomaly gradient varies with direction. Given this variation, a natural question is how the magnetic anomaly gradient affects positioning accuracy as a function of heading. This effect can be determined by first characterizing the overall positioning accuracy achievable from a given map. For example, Tkhorenko and Karshakov (2022) use Cramér-Rao Lower Bound (CRLB) analysis to evaluate how “upward continuation”of magnetic anomaly map data (i.e., stable extrapolation away from the magnetic anomaly source; (Phillips, 2005)) affect static position estimates. Similarly, Gupta et al. (2024) use a Bayesian CRLB framework to address the more general problem of directional dependency. These methods allow the magnetic anomaly map to be translated into a CRLB map, as shown in Figure 1 (right), Figure 2 (right), and Figure 3 (right). The gradient information for these maps comes from the North American Magnetic Anomaly Group (2002) and is used to calculate the Bayesian along-track CRLB maps for an aircraft traveling along a straight path at 85 m/s and 1000 m altitude at headings of 0, 60, and 120 degrees, respectively. The gradient may include errors introduced by compiling different data sources of varying quality (Saltus et al., 2023), but given that the CRLB is the lower bound on the position error, any errors in the gradient calculation are second-order effects that can be captured by the overall measurement uncertainty.

The total measurement uncertainty used to generate these CRLB plots includes the variance in the magnetometer measurement, error in the map, and calibration error introduced by electromagnetic interference from the airframe and its interaction with the core magnetic field. Of these error sources, calibration error can often dominate if sensors are not installed in areas isolated from electromagnetic interference. The calibration error can be reduced using either first-order techniques such as Tolles-Lawson calibration (Gnadt et al., 2022) or more recent machine learning-based approaches (Gnadt, 2022b; Zhai et al., 2023; Moradi et al., 2024). For example, Gnadt (2022b) suggests that calibration error can be reduced to approximately 5 nT with good sensor placement and map quality. The present study uses a more conservative calibration error of 25 nT based on recent results from the MagNav challenge data set (Gnadt et al., 2023; Nerrise et al., 2024).

The CRLB value at a location indicates the best possible navigation performance of a recursive filter, assuming of perfectly-known dynamics and a perfect inertial measurement unit providing bias-free accelerations. Because this value is necessarily a function of the gradient in a magnetic anomaly map, the expected accuracy can vary across different directions of travel. Figure 5 shows this directional sensitivity by plotting the CRLB along different headings. In this figure, CRLB values are calculated at each point using the approach presented in Gupta et al. (2024), and an ellipse is then fit through these points using least squares. Although it is possible to plot the two-dimensional CRLB covariance ellipse, the approach taken here removes any sensitivity to flight direction and therefore better represents the relationship between the direction of the magnetic anomaly gradient and the expected navigation accuracy. In this example, the ratio of principal axis standard deviations is approximately 0.4. While it may be valuable to calculate this ratio for every mapped location, such an intensive calculation is beyond the scope of this paper because it would require evaluating CRLB over a range of altitudes, groundspeeds, and headings over large maps like NAMAM (for North America) or the Earth Magnetic Anomaly Grid (EMAG2) (Maus et al., 2009) (for the entire Earth). Instead, this paper acknowledges that the positioning error statistics will certainly be mismatched and focuses on the consequences of this mismatch.

FIGURE 5

Confidence lower bound at latitude = 35.6°N, longitude = 238.1° (121.9°W)

As noted above, the standard deviations in Figure 5 are based on a measurement uncertainty of 25 nT. Even with larger uncertainty (e.g., 100 nT), the radial asymmetry shown here would persist because gradient directions are not random. Over short time intervals, if the gradient direction remains stable, the posterior variance along that direction would continue to decrease, while the variance orthogonal to the gradient would remain high.

This observation motivates the present study. Even if navigation accuracy specifications are restricted to lateral errors, MagNav can result in two deviations from the underlying assumptions:

High correlation between the longitudinal and lateral errors can cause the true radial error to exceed the given error specification, even if the lateral error meets it.
Large differences between the lateral and longitudinal error uncertainties can cause unintended separation violations in the along-track direction even when the cross-track specification is met.

This scenario is shown in Figure 6, where the shaded ellipse depicts the confidence bounds on a two-dimensional position error distribution. Error components follow a bivariate normal distribution. Even though the lateral accuracy requirement limit (dashed lines) contains the 95% cross-track uncertainty bound (dotted lines), it is still possible to violate a separation constraint with an aircraft that is slightly ahead and to the right. Conversely, calculating the probability of exceedance for a radial metric instead of a lateral or cross-track metric could inform future traffic management strategies when using MagNav.

FIGURE 6

General position error uncertainty from magnetic navigation and its effect on accuracy specifications

This work accordingly proposes a new performance metric that meets the PBN accuracy definition and compares this new metric to the current method of calculating accuracy. The paper is organized as follows. Section 2 presents the general concept of accuracy as used for current civil aviation operations. Section 3 then proposes a new metric based on the general quadratic form (Provost & Mathai, 1992) and demonstrates how this metric can be calculated. Finally, results shown in Section 4 indicate how this new metric can support tighter performance bounds under the RNP accuracy and containment region definitions.

2 BOUNDS ON ACCURACY AND THE PROBABILITY OF EXCEEDANCE

This section defines horizontal navigation accuracy in terms of the bounds on radial error when the individual error components are stochastic. To this end, let x and y denote the in-plane position errors in any coordinate frame (north and east, or along-track and cross-track). The error components x and y are typically modeled as jointly normally distributed variables with standard deviations of σ_x and σ_y respectively. For now, assume that these distributions are uncorrelated; the general problem of correlated distributions will be considered in a later section. Given these error components, the radial error r is given by $r = \sqrt{x^{2} + y^{2}}$ . This paper addresses the following question: what values of σ_x and σ_y ensure that the probability that r exceeds a threshold r_t is no greater than a specified probability of exceedance p_e?

As noted above, Section 2.3.3.3 of “ICAO Doc 9613, Performance-based Navigation (PBN) Manual” (2023) proposes two alternative methods for calculating σ_x and σ_y for a given radial distance threshold. The following two subsections describe these methods.

2.1 Radial Error Less than Threshold

The navigation accuracy requirement can be expressed in terms of the probability density function of r, a distance error threshold r_t, and the probability of exceedance p_e, as follows:

$p (r = \sqrt{x^{2} + y^{2}} > r_{t}, σ_{x}, σ_{y}) = p_{e}$ 1

In other words, the probability that the radial error r exceeds the threshold r_t should be less than or equal to p_e, with equality denoting the bounding case.

Because $\sqrt{x^{2} + y^{2}}$ is non-negative, this condition can equivalently be written as:

$p (\sqrt{x^{2} + y^{2}} \leq r_{t}, σ_{x}, σ_{y}) = 1 - p_{e}$ 2

For example, the RNP specification requires that the position error remain within the threshold for 95% of the flight time, meaning p_e = 0.05 when r_t = r_RNP.

Similarly, the containment region is defined such that the probability of exceeding r_t = r_cont is less than 10⁻⁵ or 10⁻⁷, depending on the application. This paper uses the generalized condition expressed in Equation (2) and the RNP specifications for position error.

Expressed in terms of the cumulative density function, Equation (2) becomes:

$CDF (\sqrt{x^{2} + y^{2}} = r_{t}, σ_{x}, σ_{y}) = 1 - p_{e}$ 3

where CDF(⋅) is the cumulative density function of the random variable $r = \sqrt{x^{2} + y^{2}}$ .

Assuming σ_x = σ_y, the probability density in Equation (3) corresponds to the Rayleigh distribution, which is equivalent to the chi-squared distribution with two degrees of freedom. The probability density function and CDF of the Rayleigh distribution are given by:

${PDF}_{Rayleigh} (r, σ) = \frac{r}{σ^{2}} \exp (- \frac{r^{2}}{2 σ^{2}})$ 4

${CDF}_{Rayleigh} (r, σ) = 1 - \exp (- \frac{r^{2}}{2 σ^{2}})$ 5

The general requirement for a threshold radius meeting a probability requirement under the assumptions of a Rayleigh distribution is given by:

${CDF}_{Rayleigh} (r_{t}, σ) = 1 - \exp (- \frac{r_{^{t}}^{2}}{2 σ^{2}}) = 1 - p_{e}$ 6

which can be solved for σ as shown below:

$σ = \frac{r_{t}}{\sqrt{- 2 \log p_{e}}}$ 7

Because the RNP specifications formally require that r ≤ RNP with 95% probability, Equation (7) is restated as:

$σ_{RNP} = \frac{r_{RNP}}{\sqrt{- 2 \log p_{exceedance/RNP}}} = \frac{r_{RNP}}{\sqrt{- 2 \log 0.05}} \approx 0.408 r_{RNP}$ 8

2.2 Individual Axis Errors Less than Threshold

For the second calculation method, let z denote the two-dimensional vector with components x and y. The radial error is defined as $r = \sqrt{x^{2} + y^{2}} = ‖ z ‖_{2}$ , where ∥⋅∥₂ denotes the 2-norm of the vector. According to the norm inequality between the 2-norm and the ∞-norm, ∥z∥₂≥∥z∥_∞, where ∥z∥_∞ = max(|x|, |y|). By applying a threshold condition to the individual axes, Equation (2) can be reformulated as:

$p (\max (| x |, | y |) \leq r_{t}, σ_{x}, σ_{y}) = 1 - p_{e}$ 9

Equation (9) is equivalent to:

$p (| x | \leq r_{t}, | y | \leq r_{t}, σ_{x}, σ_{y}) = 1 - p_{e}$ 10

In turn, Equation (10) allows the use of a zero-mean, uncorrelated bivariate normal distribution. with a probability density function equal to:

${PDF}_{Bivariate} (x, y, σ_{x}, σ_{y}) = \frac{1}{2 π σ_{x} σ_{y}} \exp (- \frac{x^{2}}{2 σ_{x}^{2}} - \frac{y^{2}}{2 σ_{y}^{2}})$ 11

From the distribution given in Equation (11), Equation (10) can be calculated as:

$\begin{matrix} p (| x | \leq r_{t}, | y | \leq r_{t}, σ_{x}, σ_{y}) & = & \int_{y = - r_{t}}^{y = r_{t}} \int_{x = r_{t}}^{x = r_{t}} \frac{1}{2 π σ_{x} σ_{y}} \exp (- \frac{x^{2}}{2 σ_{x}^{2}} - \frac{y^{2}}{2 σ_{y}^{2}}) d x d y \\ = & erf (\frac{r_{t}}{\sqrt{2} σ_{x}}) erf (\frac{r_{t}}{\sqrt{2} σ_{y}}) \\ = & 1 - p_{e} \end{matrix}$ 12

where erf(⋅) is the error function (Abramowitz & Stegun, 1972, p. 296).

Assuming σ_x = σ_y, Equation (12) can be solved for σ:

$σ = \frac{r_{t}}{\sqrt{2} {erf}^{-1} (\sqrt{1 - p_{e}})}$ 13

As before, for a probability of exceedance equal to 0.05, the maximum standard deviation allowed by RNP specifications is approximately:

$σ_{RNP}^{'} = \frac{r_{RNP}}{\sqrt{2} {erf}^{- 1} (\sqrt{1 - 0.05})} \approx 0.447 r_{RNP}$ 14

A comparison of Equation (14) and Equation (8) shows that applying the RNP criterion to position accuracy along individual axes allows for an approximately 9% larger standard deviation than applying it to radial accuracy. However, both approaches yield conservative bounds on the actual probability of exceedance for the radial error. This conservatism arises because these methods assume that the standard deviation of the error distribution along both axes are equal. If this assumption is violated, the larger of the two standard deviations must be used to maintain validity. As noted earlier, this assumption generally does not hold for MagNav. The next section addresses the more general problem of unequal standard deviations.

3 NAVIGATION ACCURACY USING QUADRATIC FORMS ON RANDOM VARIABLES

This section considers the general case of a two-dimensional position error vector, where x and y have a joint normal distribution. As before, the marginal distributions have standard deviations of σ_x and σ_y respectively. In this more general case, the distribution of the two-dimensional position error has a non-zero correlation coefficient ρ. The covariance matrix of the position error vector is denoted by:

$P = [\begin{matrix} σ_{x}^{2} & ρ σ_{x} σ_{y} \\ ρ σ_{x} σ_{y} & σ_{y}^{2} \end{matrix}]$ 15

Using the spectral decomposition theorem (Provost & Mathai, 1992), this covariance matrix can be written as P = VDV^⊤ where V is the matrix composed of the eigenvectors of P, and D is a diagonal matrix containing the corresponding eigenvalues. The matrix V defines a rotation from the original coordinate system in which x and y are defined to a reference frame aligned with the uncertainty ellipsoid corresponding to the normal distribution with covariance matrix P.

Let z′ = V ⁻¹z, such that z′ is the position error vector expressed in the principal axis frame. This transformed random vector is also normally distributed, with covariance given by:

$\begin{matrix} E [z^{'} {z^{'}}^{⊤}] & = & E [V^{- 1} z z^{⊤} V^{- 1 ⊤}] \\ = & V^{- 1} E [z z^{⊤}] V^{- 1 ⊤} \\ = & V^{- 1} P V^{- 1 ⊤} \\ = & V^{- 1} V D V^{⊤} V^{- 1 ⊤} \\ = & D \end{matrix}$ 16

Note that $∥ z^{'} ∥ = \sqrt{z^{' ⊤} z^{'}} = \sqrt{z^{⊤} V^{⊤} V z} = \sqrt{z^{⊤} z} = ∥ z ∥$ because V^⊤ = V⁻¹. Due to this rotation invariance, the probability density of ∥z∥ for a distribution with a positive definite covariance matrix P is equal to the probability density of ∥z′∥ for a distribution with the diagonal covariance matrix D. In other words, the variances of the new distribution are the eigenvalues of P. Without loss of generality, the distribution of ∥z∥ can therefore be considered as $‖ z ‖ = \sqrt{x^{2} + y^{2}}$ , where x and y are uncorrelated normal distributions whose standard deviations can be obtained from the eigenvalues of P (i.e. the diagonal elements of D). These eigenvalues are denoted λ_max and λ_min, respectively, and are given by:

$\begin{matrix} λ_{\max} & = & \frac{σ_{x}^{2} + σ_{y}^{2}}{2} + \frac{1}{4} \sqrt{{(σ_{x}^{2} + σ_{y}^{2})}^{2} - 4 σ_{x}^{2} σ_{y}^{2} (1 - ρ^{2})} \\ λ_{\max} & = & \frac{σ_{x}^{2} + σ_{y}^{2}}{2} - \frac{1}{4} \sqrt{{(σ_{x}^{2} + σ_{y}^{2})}^{2} - 4 σ_{x}^{2} σ_{y}^{2} (1 - ρ^{2})} \end{matrix}$ 17

Note that Equation (17) assumes that σ_x >= σ_y.

In cases where σ_x ≠ σ_y, it may be tempting to approximate the standard deviation as $σ = \sqrt{σ_{x}^{2} + σ_{y}^{2}}$ and compare this value to the threshold obtained from Equation (3). However, this approach would be incorrect because $σ_{x}^{2} + σ_{y}^{2}$ is the variance of the distribution of x + y, which can be different from the distribution of $\sqrt{x^{2} + y^{2}}$ . Most notably, x + y has a domain of (−∞, ∞), while $\sqrt{x^{2} + y^{2}}$ is restricted to [0, ∞).

In the general case where σ_x ≠ σ_y, the distribution of z follows the Hoyt distribution (Hoyt, 1947), which is also known as the Nakagami-q distribution (Tavares, 2010). (This should not be confused with the Nakagami-r distribution (Nakagami, 1960)).

The cumulative density function of the Hoyt distribution can be written as an infinite sum of incomplete gamma functions of increasing order. Ropokis et al. (2008) describe an approach for truncating this infinite sum such that the truncation error is below a predetermined threshold. In other words:

${CDF}_{Hoyt} (r) = \sum_{k = 0}^{N} \frac{c_{k}}{Γ (k + n / 2)} γ (\frac{n}{2} + k, \frac{r}{2 β})$ 18

where Γ(⋅) is the Gamma function, γ is the lower incomplete Gamma function, and n = 2 the dimensions of the system.

The scalar β is defined in Ropokis et al. (2008) as:

$β = \frac{2 \max_{i} λ_{i} \min_{i} λ_{i}}{\max_{i} λ_{i} + \min_{i} λ_{i}}$ 19

The scalar c_k is defined as follows:

$c_{k} = \frac{1}{2 k} \sum_{l = 0}^{k - 1} d_{k - l} c_{l}$ 20

where:

$c_{0} = \prod_{i = 1}^{n} {(\frac{β}{λ_{i}})}^{1 / 2}$ 21

$d_{k} = \sum_{i = 1}^{n} {(1 - \frac{β}{λ_{i}})}^{k}$ 22

For the specific case of a two-dimensional distribution, the scalar β can be expressed as:

$β = \frac{2 σ_{x}^{2} σ_{y}^{2}}{σ_{x}^{2} + σ_{y}^{2}} = \frac{2 η^{2} σ^{2}}{1 + η^{2}}$ 23

where η = σ_y/σ_x is the ratio of the smaller to the larger standard deviation (i.e., 0 ≤ η ≤ 1). Expressing the two standard deviations as a single standard deviation and scale factor more succinctly captures the relative information available along the two principal axes. When η approaches 1, the estimator provides similar uncertainty bounds along both principal axes, suggesting that the information content in the map is similar in both directions. In contrast, when η approaches 0, almost all observability is concentrated along the first principal axis.

Because calculating the inverse of the summation of a series of incomplete Gamma function is not trivial, Equation (18) is not particularly useful for deriving the standard deviation bound that meets the radial error threshold. Instead, this study uses the approximation approach developed by Tavares (2010). This method approximates the infinite sum using Gauss-Chebyshev quadrature. The number of terms needed to ensure that the approximation error does not exceed a given bound can be calculated beforehand. The Tavares approximation for the Hoyt cumulative density function is given by:

${CDF}_{Hoyt} (r) = 1 - \frac{2 q}{n (1 + q^{2})} \sum_{k = 1}^{n} \frac{\exp [- \frac{{(1 + q^{2})}^{2}}{4 q^{2}} g (q, k, n) \frac{r^{2}}{Ω}]}{g (q, k, n)} + E_{n} (q, \frac{r}{\sqrt{Ω}})$ 24

Espinosa (2019) shows that q is equal to the ratio of minimum and maximum standard deviations, as follows:

$q = \frac{\min (σ_{x}, σ_{y})}{\max (σ_{x}, σ_{y})} = η$ 25

The parameter Ω is the mean of the quadratic form (Provost & Mathai, 1992):

$Ω = σ_{x}^{2} + σ_{y}^{2} = σ^{2} (1 + η^{2})$ 26

The function g(q, k, n) is defined by Tavares (2010) as:

$\begin{matrix} g (q, k, n) & = & 1 + \frac{1 - q^{2}}{1 + q^{2}} \cos (π \frac{2 k - 1}{2 n}) \\ or, g (η, k, n) & = & 1 + \frac{1 - η^{2}}{1 + η^{2}} \cos (π \frac{2 k - 1}{2 n}) \end{matrix}$ 27

Finally, E_n is the error associated with truncating the series at n terms, given by:

$| E_{n} (q, \frac{r}{\sqrt{Ω}}) | \leq | E_{n} (q, 0) | = \frac{2}{{(\frac{1 + q}{1 - q})}^{2 n} + 1} = \frac{2}{{(\frac{1 + η}{1 - η})}^{2 n} + 1}$ 28

From Equation (25), η is known beforehand (typically from the estimator), so the number of terms n required to meet a given tolerance ε can be obtained by solving Equation (28):

$n = ceil [\frac{1}{2} \log (\frac{2}{ε} - 1) / \log (\frac{1 + η}{1 - η})]$ 29

A reasonable choice for the tolerance ε is some fraction of the exceedance probability p_e. Once n has been determined from Equation (29), the following expression must be solved to determine the standard deviation σ that satisfies a given containment radius and exceedance probability:

${CDF}_{Hoyt} (r_{t}, σ, η) = 1 - p_{e}$ 30

This equation can be solved for σ using the truncated sum of terms in Equation (24):

$\begin{matrix} \frac{2 η}{n (1 + η^{2})} \sum_{k = 1}^{n} \frac{\exp [- \frac{{(1 + η^{2})}^{2}}{4 η^{2}} g (η, k, n) w]}{g (η, k, n)} = p_{e} \\ w = \frac{r_{t}^{2}}{σ^{2} (1 + η^{2})} \end{matrix}$ 31

Because Equation (31) is nonlinear in σ, it is solved using Newton-Raphson iteration:

Initialize $\hat{w}$ with the variance calculated from the Rayleigh threshold in Equation (7):
$\hat{w} = \frac{r_{t}^{2}}{σ^{2} (1 + η^{2})} = \frac{2 \log p_{e}}{1 + η^{2}}$ 32
Calculate the function $f (\hat{w})$ given by Equation (31):
$f (\hat{w}) = \frac{2 η}{n (1 + η^{2})} \sum_{k = 1}^{n} \frac{\exp [- \frac{{(1 + η^{2})}^{2}}{4 η^{2}} g (η . k, n) \hat{w}]}{g (η . k, n)}$ 33
Calculate $f^{'} (\hat{w})$ , which is the derivative of f(w) with respect to w evaluated at $\hat{w}$ :
$f^{'} (\hat{w}) = \frac{\partial f (w)}{\partial w} |_{w = \hat{w}} = - \frac{(1 + η^{2})}{2 n η} \sum_{k = 1}^{n} \exp [- \frac{(1 + η^{2})}{4 η^{2}} g (η, k, n) \hat{w}]$ 34
Calculate a correction factor ∆w:
$Δ w = \frac{p_{e} - f (\hat{w})}{f^{'} (\hat{w})}$ 35
Update the measurement: ${\hat{w}}^{'} = \hat{w} + Δ w$
If $| f ({\hat{w}}^{'}) - f (\hat{w}) | > ε$ , update $\hat{w} \leftarrow {\hat{w}}^{'}$ and repeat Steps (2) through (5). When $| f ({\hat{w}}^{'}) - f (\hat{w}) | \leq ε$ , calculate σ from Equation (31):
$σ = \sqrt{\frac{r_{t}^{2}}{\hat{w} (1 + η^{2})}}$ 36

This method of computing uncertainty bounds from the Hoyt distribution is generally stable except for limiting cases when η → 1 or η → 0. These two limiting cases are discussed in the following subsections.

3.1 Limiting Case: η → 1

When η = 1, the denominator in Equation (29) becomes indeterminate; however, the limit still exists and yields n = 1. With n = 1 and q = η = 1, the Gauss-Chebyshev quadrature form of Equation (24) is identical to the Rayleigh distribution in Equation (5). Therefore Equation (7) can be used for this specific limiting case.

For a given tolerance ε, the number of terms n is a decreasing function of η with a minimum value of 1. The calculation of the number of terms in Gauss-Chebyshev quadrature can be safeguarded by ensuring that n is at least equal to 1 when η is greater than some threshold value η_max. This threshold can be determined by solving Equation (29) with n = 1:

$\begin{matrix} \log (\frac{1 + η_{\max}}{1 - η_{\max}}) = \frac{1}{2} \log (\frac{2}{ε} - 1) \\ or, & \frac{1 + η_{\max}}{1 - η_{\max}} = \exp [\frac{1}{2} \log (\frac{2}{ε} - 1)] \\ \Rightarrow & η_{\max} = \frac{\exp [\frac{1}{2} \log (\frac{2}{ε} - 1)] - 1}{exp [\frac{1}{2} \log (\frac{2}{ε} - 1)] + 1} = \frac{\sqrt{\frac{2}{ε} - 1 - 1}}{\sqrt{\frac{2}{ε} - 1 + 1}} \end{matrix}$ 37

For example, when the required tolerance is ε = 10⁻⁹, Equation (37) gives a limiting ratio of η_max = 0.9999553. If η is larger than this value, it is appropriate to use the Rayleigh distribution without concern for overestimating the bound within the specified accuracy limits.

3.2 Limiting Case: η → 0

As η approaches 0, the number of terms required to calculate the sum in Equation (24) approaches infinity because the denominator in Equation (29) approaches zero. In the limiting case where η = 0, the two-dimensional distribution reduces to a one-dimensional distribution where the threshold requirement is given by CDF(|x| ≤ r_t) = 1 − p_e. This cumulative density function corresponds to the one-sided (or half-normal) Gaussian distribution and can be derived from the cumulative density function of the standard two-sided Gaussian distribution. In fact, the value of σ can be calculated by letting σ_y → ∞ in Equation (12), for which the error function approaches unity. Consequently:

$e r f (\frac{r_{t}}{\sqrt{2 σ}}) = 1 - p_{e}$ 38

which can be solved as follows:

$σ = \frac{r_{t}}{\sqrt{2} e r f^{- 1} (1 - p_{e})}$ 39

Let n_max be the maximum number of terms allowed in the quadrature formula. For a given error tolerance ε, the lower bound on η, denoted by η_min can be obtained by solving Equation (29). This value is given by:

$η_{\min} = \frac{(\frac{2}{ε} - 1)^{1 / (2 n_{\max})} - 1}{(\frac{2}{ε} - 1)^{1 / (2 n_{\max})} + 1}$ 40

Therefore, if η < η_min, the one-sided Gaussian solution in Equation (39) should be used.

4 RESULTS AND ANALYSIS

Figure 7 compares the one-sided Gaussian (circles), Hoyt (squares), and Rayleigh distributions (diamonds) by plotting the probability of exceedance (i.e., the survival function, or one minus the cumulative density) on a logarithmic scale. As discussed earlier, the one-sided Gaussian is the limiting distribution when the ratio of standard deviations along the minor and major axes approaches zero. Conversely, the Rayleigh distribution is the limiting distribution when the ratio of standard deviations approaches one (i.e. both standard deviations are equal). For Figure 7, η = 0.8. The figure also shows that the Hoyt distribution calculated using quadratures agrees closely with the version calculated using the truncated sum of weighted incomplete gamma functions. The difference between the two is on the order of 10⁻¹³, as depicted by the dashed line.

FIGURE 7

Comparison of probability of exceedance across the one-sided Gaussian, Hoyt, and Rayleigh distributions

Table 1 presents the thresholding standard deviation σ_{limit_Hoyt} calculated from the Hoyt distribution for r_t = r_RNP, corresponding to a probability of exceedance of 0.05. The implementation outlined in Equations (24) through (36) is used to calculate this value to an error tolerance of 10⁻¹². Standard deviations shown in Table 1 are scaled by r_RNP to nondimensionalize the result. In the third column of Table 1, the limiting standard deviation from the Hoyt distribution is divided by the standard deviation of the Rayleigh distribution, which is equal to 0.408 r_RNP as shown in Equation (8). This ratio demonstrates the potential benefit of using a Hoyt distribution for radial error characterization in MagNav. For example, if the standard deviation along the minor axis of the covariance ellipse is an order of magnitude smaller than the standard deviation along the major axis (i.e., η ≈ 0.1), the benefit becomes significant. In general, even when the principal axis standard deviations differ only moderately (e.g. the ratio of the smaller to the larger standard deviation is 0.5 or lower), using the radial error distribution instead of a single-axis error distribution allows for a 20% relaxation in the required standard deviation for a 95th percentile bound. This added flexibility can support navigation filter designs by removing the need for precise tuning on each axis and instead requiring tuning to reduce the radial error only. Conversely, applying a 95% accuracy requirement to the radial error leads to a stronger guarantee on separation assurance than would be achieved with single axis-based constraints.

View this table:

TABLE 1

RNP Requirement Relaxation with Hoyt Distribution.

An alternative approach involves using a modified version of Equation (12) where σ_x = σ and σ_y = ησ:

$e r f (\frac{r_{t}}{\sqrt{2 σ}}) e r f (\frac{r_{t}}{\sqrt{2 η σ}}) = 1 - p_{e}$ 41

This equation can be solved for σ via Newton-Raphson iteration, but it assumes that numerical implementations of the error function and its inverse are available. However, Equation (41) applies a constraint on the infinity norm of the two-dimensional position error, not on the radial error (i.e., 2-norm of the error). For a two-dimensional random variable z, the region defined by ∥z∥_∞ ≤ r_t is always greater than the region defined by ∥z∥₂ ≤ r_t. Consequently, p(∥z∥_∞ ≤ r_t) ≥ p(∥z∥₂ ≤ r_t), and the probability of exceeding r_t follows p(∥z∥_∞ > r_t) < p(∥z∥₂ > r_t). This means that enforcing the RNP accuracy condition on the infinity norm does not provide any meaningful bound on the radial error.

Figure 8 generalizes the results from Table 1 to a range of exceedance probabilities. In this figure, level curves are plotted with the x-axis representing the ratio of principal major and minor axis stand deviations and the y-axis showing the log-transformed probability of exceedance. These curves show that, as η increases, the bounding standard deviation from the Hoyt distribution approaches the value obtained from the Rayleigh distribution, as expected. This convergence means that, for sufficiently high η, the penalty incurred by using a Rayleigh approximation is negligible, making it an acceptable alternative to the Hoyt distribution. The same holds for sufficiently low probabilities of exceedance. For example, containment region statistics that require a probability of exceedance of 10⁻⁷ can still use a Rayleigh approximation, which only overestimates the limiting σ value by 5%.

FIGURE 8

Ratio of standard deviation limits obtained from Hoyt and Rayleigh distributions for different probabilities of exceedance

Finally, Figure 9 shows the number of Newton-Raphson iterations required to calculate the bounding standard deviation when the required tolerance is 10⁻¹². The required number of iterations is a function of both η and the target exceedance probability, but in all cases, convergence occurs within 12 iterations. This number is reduced to 8 if the required tolerance is reduced to 10⁻¹⁰. If computational limitations constrain the number of possible iterations, the gradient calculation in Equation (34) can be used to obtain a first-order correction to the Rayleigh distribution. However, this paper does not evaluate the approximation error due to this limit.

FIGURE 9

Required number of Newton-Raphson iterations

5 CONCLUSIONS

Integrating MagNav into civil aviation operations requires adapting concepts originally developed for GPS-basd navigation, such as RNP. The definitions of RNP and related accuracy metrics are rooted in statistical analyses where the probability of exceedance is derived from the standard deviations of one-dimensional normal distributions. However, MagNav does not always satisfy the underlying assumptions used to define these accuracy metrics. To address this problem, this paper quantified the errors incurred while making these assumptions. The paper then developed a radial error model, based on the Hoyt distribution, that more rigorously adheres to the spirit of RNP and similar navigation accuracy definitions. Because the probability density function, cumulative density function, and inverse functions of the Hoyt distribution are not easily expressed using elementary functions, the paper also developed robust numerical methods to evaluate these functions and make their computation tractable. Given the robust, explainable, and computationally tractable nature of this approach, the radial error statistic presented here is well-suited for quantifying navigation accuracy in civil aviation.

HOW TO CITE THIS ARTICLE

Sengupta, P. (2025). A horizontal accuracy metric for magnetic navigation. NAVIGATION, 72(4). https://doi.org/10.33012/navi.717

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

↵
Abramowitz, M., & Stegun, I. A. (1972). Handbook of mathematical functions with formulas, graphs, and mathematical tables (10th ed.). Marcel Dekker. https://dl.acm.org/doi/10.5555/1098650
Google Scholar
↵
Blanch, J., Walter, T., Enge, P., Lee, Y., Pervan, B., Rippl, M., Spletter, A., & Kropp, V. (2015). Baseline advanced RAIM user algorithm and possible improvements. IEEE Transactions on Aerospace and Electronic Systems, 51(1), 713–732. https://doi.org/10.1109/TAES.2014.130739
Google Scholar
↵
Canciani, A. J. (2022). Magnetic navigation on an F-16 aircraft using online calibration. IEEE Transactions on Aerospace and Electronic Systems, 58(1), 420–434. https://doi.org/10.1109/TAES.2021.3101567
Google Scholar
↵
Canciani, A. J., & Raquet, J. F. (2016). Absolute positioning using the Earth’s magnetic anomaly field. NAVIGATION, 63(2), 111–126. https://doi.org/10.1002/navi.138
Google Scholar
↵
Crassidis, J. L., & Junkins, J. L. (2004). Optimal estimation of dynamical systems. Chapman & Hall.
Google Scholar
↵
Espinosa, P. R. (2019). Analysis of Gaussian quadratic forms with application to statistical channel modeling [PhD Thesis, Universidad de Málaga]. https://hdl.handle.net/10630/19450
Google Scholar
↵
Federal Aviation Administration. (2024). 2024 Federal aviation regulations/aeronautical information manual. Aviation Supplies & Academics. https://www.faa.gov/air_traffic/publications/atpubs/aim_html/
Google Scholar
↵
Gnadt, A. R. (2022a). Advanced aeromagnetic compensation models for airborne magnetic anomaly navigation [PhD Thesis, Massachusetts Institute of Technology]. https://hdl.handle.net/1721.1/145137
Google Scholar
↵
Gnadt, A. R. (2022b). Machine learning-enhanced magnetic calibration for airborne magnetic anomaly navigation. Proc. of the AIAA SCITECH 2022 Forum, San Diego, CA. https://doi.org/10.2514/6.2022-1760
Google Scholar
↵
Gnadt, A. R., Belarge, J., Canciani, A., Carl, G., Conger, L., Curro, J., Edelman, A., Morales, P., Nielsen, A. P., O’Keeffe, M. F., Rackauckas, C. V., Taylor, J., & Wollaber, A. B. (2023). Signal enhancement for magnetic navigation challenge problem. https://arxiv.org/abs/2007.12158
Google Scholar
Gnadt, A. R., Wollaber, A. B., & Nielsen, A. P. (2022). Derivation and extensions of the Tolles-Lawson model for aeromagnetic compensation. https://arxiv.org/abs/2212.09899
Google Scholar
↵
Gupta, A., Sengupta, P., Phernetton, R., & Sosanya, A. (2024). Lower bounds on magnetic navigation performance as a function of magnetic anomaly map quality. Proc. of the 2024 AIAA DATC/IEEE 43rd Digital Avionics Systems Conference (DASC), San Diego, CA, 1–7. https://doi.org/10.1109/DASC62030.2024.10749342
Google Scholar
↵
Hoyt, R. S. (1947). Probability functions for the modulus and angle of the normal complex variate. The Bell System Technical Journal, 26(2), 318–359. https://doi.org/10.1002/j.1538-7305.1947.tb01318.x
Google Scholar
↵
International Civil Aviation Organization. (2023). ICAO doc 9613, performance-based navigation (PBN) manual (5th ed.). https://pbnportal.eu/dam/jcr:ca055ef7-5fa7-45e1-b4ee-7319dbc486b6/9613_unedited_en_V5.pdf
Google Scholar
↵
Langley, R. B. (1999). Dilution of precision. GPS World, 10(5), 52–59. https://gge.ext.unb.ca/Resources/gpsworld.may99.pdf
Google Scholar
↵
Maus, S., Barckhausen, U., Berkenbosch, H., Bournas, N., Brozena, J., Childers, V., Dostaler, F., Fairhead, J. D., Finn, C., von Frese, R. R. B., Gaina, C., Golynsky, S., Kucks, R., Lühr, H., Milligan, P., Mogren, S., Müller, R. D., Olesen, O., Pilkington, M., … Caratori Tontini, F. (2009). EMAG2: A 2–arc min resolution Earth magnetic anomaly grid compiled from satellite, airborne, and marine magnetic measurements. Geochemistry, Geophysics, Geosystems, 10(8). https://doi.org/10.1029/2009GC002471
Google Scholar
↵
Moradi, M., Zhai, Z.-M., Nielsen, A., & Lai, Y.-C. (2024). Random forests for detecting weak signals and extracting physical information: A case study of magnetic navigation. APL Machine Learning, 2(1), 016118. https://doi.org/10.1063/5.0189564
Google Scholar
↵
Nakagami, M. (1960). The m-distribution — A general formula of intensity distribution of rapid fading. In W. C. Hoffman (Ed.), Statistical methods in radio wave propagation, 3–36. Pergamon. https://doi.org/10.1016/B978-0-08-009306-2.50005-4
Google Scholar
↵
Nerrise, F., Sosanya, A. S., & Neary, P. (2024). Physics-informed calibration of aeromagnetic compensation in magnetic navigation systems using liquid time-constant networks. Proc. of the 37th Annual Conference on Neural Information Processing Systems (NeurIPS 2023), New Orleans, LA. https://neurips.cc/media/neurips-2023/Slides/76255.pdf
Google Scholar
↵
North American Magnetic Anomaly Group. (2002). Magnetic anomaly map of North America. U.S. Geological Survey. https://doi.org/10.3133/70211067
Google Scholar
↵
Phillips, J. D. (2005). Potential-field continuation: Past practice vs. modern methods. In SEG technical program expanded abstracts 1996, 1411–1414. https://doi.org/10.1190/1.1826376
Google Scholar
↵
Provost, S. B., & Mathai, A. M. (1992). Quadratic forms in random variables: Theory and applications. Marcel Dekker.
Google Scholar
↵
Ropokis, G. A., Rontogiannis, A. A., & Mathiopoulos, P. T. (2008). Quadratic forms in normal RVs: Theory and applications to OSTBC over Hoyt fading channels. IEEE Transactions on Wireless Communications, 7(12), 5009–5019. https://doi.org/10.1109/T-WC.2008.070830
Google Scholar
↵
Saltus, R. W., Chulliat, A., Meyer, B., Bates, M., & Sirohey, A. (2023). Magnetic anomaly grid and associated uncertainty from marine trackline data: The Caribbean alternative navigation reference experiment (CANREx). Earth and Space Science, 10(11). https://doi.org/10.1029/2023EA002958
Google Scholar
↵
SC-159. (2020a). DO229F: Minimum operational performance standards (MOPS) for global positioning system/satellite-based augmentation system airborne equipment. RTCA, Inc. https://my.rtca.org/productdetails?id=a1B1R0000092ubbUAA
Google Scholar
↵
SC-159. (2020b). DO384: Minimum operational performance standards (MOPS) for GNSS aided inertial systems. RTCA, Inc. https://my.rtca.org/productdetails?id=a1B1R00000LowzIUAR
Google Scholar
↵
Tavares, G. N. (2010). Efficient computation of Hoyt cumulative distribution function. Electronics Letters, 46(7), 537–539. https://doi.org/10.1049/el.2010.0189
Google Scholar
↵
Tkhorenko, M., & Karshakov, E. (2022). Estimating the potential accuracy of magnetic navigation based on magnetic survey data. Proc. of the 2022 29th Saint Petersburg International Conference on Integrated Navigation Systems (ICINS), Saint Petersburg, Russian Federation, 1–3. https://doi.org/10.23919/ICINS51784.2022.9815357
Google Scholar
↵
Zhai, Z.-M., Moradi, M., Kong, L.-W., & Lai, Y.-C. (2023). Detecting weak physical signal from noise: A machine-learning approach with applications to magnetic-anomaly-guided navigation. Physical Review Applied, 19, 034030. https://doi.org/10.1103/PhysRevApplied.19.034030
Google Scholar

Previous Next

A Horizontal Accuracy Metric for Magnetic Navigation

Abstract

1 INTRODUCTION