Cramér–Rao Lower Bound for Magnetic Field Localization around Elementary Structures

The determination of a mobile terminal’s position with high accuracy and ubiquitous coverage is still challenging. Global satellite navigation systems (GNSSs) provide sufficient accuracy in areas with a clear view to the sky. For GNSS-denied environments like indoors, complementary positioning technologies are required. A promising approach is to use the Earth’s magnetic field for positioning. In open areas, the Earth’s magnetic field is almost homogeneous, which makes it possible to determine the orientation of a mobile device using a compass. In more complex environments like indoors, ferromagnetic materials cause distortions of the Earth’s magnetic field. A compass usually fails in such areas. However, these magnetic distortions are location dependent and therefore can be used for positioning. In this paper, we investigate the influence of elementary structures, in particular a sphere and a cylinder, on the achievable accuracy of magnetic positioning methods. In a first step, we analytically calculate the magnetic field around a sphere and a cylinder in an outer homogeneous magnetic field. Assuming a noisy magnetic field sensor, we investigate the achievable positioning accuracy when observing these resulting fields. For our analysis, we calculate the Cramér–Rao lower bound, which is a fundamental lower bound on the variance of an unbiased estimator. The results of our investigations show the dependency of the positioning error variance on the magnetic sensor properties, in particular the sensor noise variance and the material properties, i.e., the relative permeability of the sphere with respect to the cylinder and the location of the sensor relative to the sphere with respect to the cylinder. The insights provided in this work make it possible to evaluate experimental results from a theoretical perspective.


Introduction
Position information is becoming increasingly important in mobile communication systems.Location-based services and navigation applications require such position information with sufficient accuracy and availability in both outdoor and indoor environments.In addition to applications, lower layers in mobile communication systems can benefit from information about the mobile terminals' locations as well [1].
There are a variety of technologies for the determination of a mobile terminal's position.Global navigation satellite systems [2] like GPS, GLONASS, Galileo, Beidou, etc. provide good accuracy and coverage in outdoor areas.Indoors, however, non-line-ofsight and multipath propagation, together with weak signal power levels at the receiver, make satellite-based positioning challenging or even impossible.In such challenging environments, positioning based on terrestrial wireless networks is considered to be the complementary technology.Received signals in a terrestrial wireless network provide a relation between the position of a mobile terminal and the positions of the base stations in the network.Several characteristics of a received signal carry information about the position of a mobile terminal [3].The time of arrival (TOA) and time difference of arrival (TDOA) exploit the signal propagation delay between base stations and the mobile terminal for positioning.Other methods make use of the angle of arrival (AOA) or the received signal strength (RSS) for mobile terminal positioning.Compared to satellite signals, terrestrial mobile radio signals provide higher power levels at the receiver.As in the case of satellite-based positioning, multipath and non-line-of-sight propagation also degrade the positioning performance based on terrestrial wireless networks severely.Furthermore, at least three base stations must be observable to a mobile terminal in order to obtain a two-dimensional position fix.This is hardly achieved in today's mobile radio networks.In the future, this situation might either change due to denser network cell structures, or the requirement of receiving multiple base stations will become obsolete by applying cooperative positioning methods [4].
Especially for pedestrian indoor navigation, the use of inertial sensors became popular in recent years.Such sensors are cheap and widely deployed in today's mobile devices.However, if conventional dead reckoning is applied using inertial sensors, the positioning error grows cubically in time due to sensor drift.In [5], the author proposes a pedestrian tracking system using a foot-mounted inertial sensor.By detecting the resting phase of the foot during a typical human walk, so-called zero-velocity updates are applied, which reduce the positioning error from cubic to linear.However, the problem of an unlimited growth of the positioning error remains.A step forward for solving the drift problem is to use landmarks, i.e., location-dependent characteristics that can be observed and identified.By regularly revisiting such landmarks, remaining inertial sensor drifts can be compensated.Mobile radio networks or dedicated positioning systems can provide such landmarks.To achieve an appropriate local distribution of landmarks, i.e., reference or base stations, some installation of additional infrastructure might be required.Additionally, this infrastructure has to be calibrated so that the positions of these reference stations are known.Both the installation and calibration of positioning infrastructure is time-consuming and costly.Therefore, a common approach is to use landmarks that can be detected with today's mobile terminal low-cost sensor equipment.Using the principles of simultaneous localization and mapping (SLAM), the positions of the mobile terminal and the landmarks are estimated at the same time, i.e., calibration happens on the fly.In [6], the authors propose a pedestrian navigation system that combines odometry, obtained from a foot-mounted inertial sensor, and RSS measurements of WiFi base stations as landmarks.Both the mobile terminal and the WiFi bass stations' positions are estimated using a SLAM approach based on Bayesian filtering.The authors of [7] introduce an indoor positioning principle, solely based on a foot mounted inertial sensor, called FootSLAM.The structure of walking routes are used as location-dependent characteristics.The FootSLAM algorithm method detects and maps these walking routes.When revisited, mapped walking routes are recognized that are used for the compensation of the inertial sensor drift.
Another kind of location-dependent characteristic that can be used for positioning is the Earth's magnetic field.Ferromagnetic materials placed in that field cause distortions.Regardless if being distorted or not, the Earth's magnetic field observed by a magnetic sensor in a mobile device in general depends on the position and/or the orientation of the mobile device.A variety of positioning methods based on the observation of the Earth's magnetic field have been proposed and studied, which relate to our investigations.Many of these methods come with experimental validation.This paper aims to provide basic theoretical dependencies for magnetic field-based positioning so that the experimental results regarding positioning performance can ultimately be evaluated and explained from a theoretical perspective.

Related Work
The discovery and observation of magnetic phenomena goes back even before the common era.As one result, the compass had already been developed centuries ago [8].A compass indicates the direction of a surrounding magnetic field.In open areas, the Earth's magnetic field provides orientation towards the magnetic poles of the Earth.More complex environments like indoors show severe distortions of the Earth's magnetic field, which are mainly caused by surrounding ferromagnetic materials such as reinforcing steel, for instance.A compass usually fails in such areas.Such distortions, however, are location-dependent characteristics and can be considered as a local magnetic fingerprint.
In [9], the authors study the feasibility of using such fingerprints for positioning.They found that the Earth's magnetic field is temporarily stable enough, and its distortions can be sufficiently significant in order to achieve submeter or even down to decimeter positioning accuracy with fingerprinting methods.Obviously, higher local distortions in the magnetic field provide higher positioning accuracy.Fingerprinting positioning methods require a database, which must be created in a training or calibration phase.
Investigations in [10] address the behavior of the Earth's magnetic field in the vicinity of ferromagnetic materials.It was shown that ferromagnetic objects like pillars can be identified and used as landmarks by observing significant changes in the magnetic field strength around them.Furthermore, the experiments described in [10] show a good match between the measured magnetic field strength and its prediction by analytical approximations.Also, long-term stability and reproducibility of the magnetic field observations are found to be suitable for magnetic fingerprinting positioning.
In [11], Angermann et al. proposed to exploit the characteristic Earth magnetic field perturbations for indoor positioning.The authors discuss various methods for the measurement and mapping of the Earth's magnetic field with sufficient granularity.For repeated magnetic field measurements, only little noise was identified, which indicated good stability and reproducibility of the Earth magnetic field observations.An autonomous robotic platform for high resolution and complete mapping of the Earth's magnetic field has been developed and described in [12].
The suitability of the magnetic field was also shown for road [13] and railway environments [14].Even in the airspace, the magnetic field shows distortions that, in combination with an inertial navigation system, can be exploited for position estimation [15].Unlike in the other environments, the distortions in the airspace are not caused by nearby magnetic material but by the crustal field of the Earth [15], which typically is smaller than the distortions observed, e.g., in indoor environments.
In [16], the theoretically achievable accuracy of magnetic localization for a wheeled robot driving through an indoor environment has been analyzed by deriving a Bayesian lower bound.The bound was based on a Gaussian process fitted to a data set of densely measured magnetic field vectors.Thus, the bound requires an exhaustive measurement campaign to be calculated, but once the data set is available, the bound can be used to assess the quality of different filter algorithms.The paper showed that particle filters perform well but are not optimal with respect to to their mean square error.
In contrast to performing magnetic field mapping as a required calibration step prior to solving the actual positioning task, several research works propose using SLAM.In [17], the authors provide an experimental proof of a concept for a SLAM approach that uses the Earth's magnetic field anomalies in combination with the odometry data of a robotic platform.The magnetic field there was parametrically modeled using Gaussian processes.The combination of magnetic and inertial sensors in a SLAM-based positioning solution has been proposed in [18][19][20].With this sensor combination, the authors of [19] obtained positioning errors in the order of 10 cm to 20 cm when walking through a large building.
The works and methods mentioned above utilize the distortions of the Earth's magnetic field for positioning.They mainly rely on experimental work for characterizing and mapping the Earth's magnetic field and its perturbations for positioning-either with a traditional calibration step prior to running positioning algorithms or with SLAM-based solutions.In order to quantify the positioning performance, estimated positions obtained during experiments have been compared to the corresponding ground truth.Our work presented in this paper complements these results by analyzing the achievable positional accuracy using an estimation-theoretic approach based on the computation of Fisher information and the Cramér-Rao lower bound.

Contribution
In this paper, we aim to lower bound the achievable error variance for positioning methods based on the observation of the ambient magnetic field.In particular, we address the following questions: • How does the achievable positioning performance depend on the magnetic sensor properties, i.e., the sensor noise variance?• How does the positioning performance depend on the distance from ferromagnetic objects that cause magnetic field distortions?• How does the ferromagnetic material properties, in particular the relative permeability, influence the achievable positioning performance?
To answer these questions, we provide an analytical description of the magnetic field around a sphere and a cylinder placed in an outer homogeneous static magnetic field in Sections 2 and 3. We use these analytical field descriptions and calculate the Cramér-Rao lower bound for position estimation based on noisy magnetic field observations in Section 4. Section 5 provides examples and evaluations of the achievable positioning performance for a ferromagnetic cylinder in the Earth's magnetic field.

The Static Magnetic Field
Maxwell's Equations-a set of four partial differential equations-describe the behavior and interrelation of charges, currents, and electric and magnetic fields.According to these equations, dynamic electric and magnetic fields influence each other through nonzero partial derivatives of those fields with respect to time.In the static case, however, these derivatives are zero, and Maxwell's system of equations decouples.So, we obtain two equations-in particular the Maxwell-Ampère equation and Gauss's law for magnetism-which describe any static magnetic field.

Maxwell's Equations for the Static Magnetic Field
Maxwell's equations for the static magnetic field can be found in many textbooks on electromagnetic field theory, e.g., in [21,22].The Maxwell-Ampère equation for the static magnetic strength vector field H in the absence of currents reads as rot H = 0. ( Another equation, describing the magnetic field, is Gauss's law for magnetism.It states that the magnetic flux density vector field is defined as with B = µ H being is a divergence-free vector field.In other words, this means that there are no magnetic charges.Here, µ is a constant describing the material's permeability.With these conditions, we can express the magnetic strength vector field as the gradient of a scalar potential field Ψ.This choice fulfills (1), since the curl of a gradient field is zero in general.Substituting (3) into (2) yields Laplace's equation: This partial differential equation usually comes with boundary conditions.Thus, a solution to that boundary value problem solves Laplace's equation and additionally satisfies the boundary conditions.

Boundary Values and Conditions
As mentioned above, we have to find a scalar magnetic potential field that is a solution of Laplace's equation and satisfies given boundary values and conditions.General bound-ary value constraints are the continuity of the parallel magnetic strength vector components H ∥ i and the perpendicular magnetic flux density vector components B ⊥ i at the boundary surfaces between materials, indexed by i, with different permeabilities, i.e., The boundary condition for the magnetic field strength components (5) follow from the Maxwell-Ampère Equation (1).The continuity of the perpendicular magnetic flux density ( 6) is a consequence of Gauss's law for magnetism (2).
Subsequently, we investigate the magnetic field that results from the presence of a ferromagnetic body like an iron sphere or cylinder in an outer magnetic field.It is reasonable to assume that the influence of the ferromagnetic body on the magnetic field decays with increasing distance.This yields the boundary value at infinite distance to the ferromagnetic body in such a way that the solution of the boundary value problem, −∇Ψ, converges to the original field H ∞ , which would be observable in absence of the ferromagnetic structure.This assumption leads to the Neumann boundary condition: Equivalently, we can require the magnetic potential field to take on specific values.In this case, we obtain the Dirichlet boundary condition which, except for an additive constant, is equivalent to the Neumann boundary condition.
The usual procedure of solution for boundary value problems is to find a general solution for Laplace's equation for each volume having a different permeability.Using their remaining degrees of freedom, these solutions are then aligned to each other in order to meet the given boundary values and constraints.

Ferromagnetic Structures in a Homogeneous Magnetic Field
As mentioned earlier, ferromagnetic materials cause distortions in the Earth's ambient magnetic field, which are inherently location dependent and can thus be used for position estimation.In daily life, there are many environments, like indoors, where we find ferromagnetic materials in close proximity.Just think of the reinforcement steel in the concrete hulls of buildings as an example.
In the following sections, we provide an analytical description of the magnetic potential field or magnetic field strength for a sphere and a cylinder, respectively, which we place into an external homogeneous magnetic field.

Sphere
We consider an outer homogeneous magnetic strength vector field with magnitude H ∞ , which is oriented in z direction of a Cartesian coordinate system.The second part of (9) provides the description of this field in spherical coordinates with unit vectors e r in the radial direction and e θ in the polar direction.This static magnetic field can be obtained from the corresponding scalar magnetic potential field: Now, let us place a sphere with a permeability µ in and a radius r sph into this magnetic field at the origin of the coordinate system, as shown in Figure 1.The resulting magnetic field can be calculated by solving Laplace's Equation ( 4) for the magnetic scalar potential field.We split this field into the magnetic potential field Ψ in within the sphere and the outer field Ψ out .At the sphere's surface, the boundary conditions ( 5) and ( 6) must hold.Additionally, the outer field has to fulfill the boundary condition lim r→∞ Ψ out = Ψ ∞ .
For solving Laplace boundary value problems, a common approach is to apply the 'separation of variables' method.
There, we assume a factorized solution of Ψ(r, φ, θ) = R(r) Φ(φ) Θ(θ), where each of the three factors only depends on one coordinate.This approach separates the Laplace differential equation into three ordinary differential equations-one for each factor.These ordinary differential equations are easier to solve.The solution of the Laplace boundary value problem through the separation of variables is well known and therefore not in the focus of this paper.Interested readers are referred to textbooks, e.g., ([21], p. 247 f.) or Appendix B. At the end, we obtain a solution of the Laplace equation for the outer magnetic potential field of our problem as Correspondingly, the inner magnetic potential field is With the gradient operator ∂ ∂φ for spherical coordinates, we obtain the corresponding magnetic field strengths as and For convenience, we have set µ r = µ in µ out .It is straightforward to verify this solution by insertion into (4)- (8).For further investigations, we will require the magnetic strength vector field outside the sphere in Cartesian coordinates.From (13) and with e r = sin(θ) cos(φ) e x + sin(θ) sin(φ) e y + cos(θ) e z (15) we obtain Figure 2 shows the magnetic field strength for a ferromagnetic sphere with a relative permeability of µ r = µ in µ out = 100 as a color plot exemplary for slice planes y = 0 and z = 0.The magnetic field strength is expressed by its magnitude normalized to the homogeneous field strength H ∞ .Its orientation is shown by magnetic flux lines.

Cylinder
We assume a cylinder of infinite length with a permeability µ in and a radius r sph located along the z axis of a cylindrical coordinate system, as shown in Figure 3. Due to its permeability µ in being different from that of the outer medium µ out , this cylinder influences an outer homogeneous magnetic field.First, let us consider a homogeneous outer magnetic field with arbitrary direction.Since the dimension of the cylinder in the z direction is infinite, there is no dependency of the resulting magnetic field on dimension z.For this reason, it is obvious that the cylinder does not influence the z component of an outer homogeneous magnetic field.Therefore, we consider a homogeneous outer magnetic field strength with its direction in the x-y plane.Without loss of generality, we assume a homogeneous magnetic strength vector field in the x direction with magnitude H ∞ , which we can express as The resulting magnetic field can be obtained by solving Laplace's Equation ( 4) for the magnetic scalar potential field.We split this field into the magnetic potential field Ψ in within the cylinder and the outer field Ψ out .The outer field has to fulfill the boundary condition − grad Ψ out = H ∞ e x , r → ∞.The solution for the outer field is Correspondingly, the inner magnetic potential field is A derivation for this solution can be found in Appendix C. With the gradient operator ∂z for cylindrical coordinates, we obtain the corresponding magnetic field strengths as outside the cylinder and inside the cylinder.For convenience, we have set µ r = µ in µ out .For further investigations, we express the magnetic strength vector field in Cartesian coordinates.From (21) and with e r = cos(φ) e x + sin(φ) e y (23) we obtain Figure 4 shows the magnetic field strength for a ferromagnetic sphere with a relative permeability of µ r = µ in µ out = 100.The magnetic field strength is expressed by its magnitude normalized to the homogeneous field strength H ∞ .Its orientation is shown by magnetic flux lines.Since the magnetic field components neither contain components in the z direction nor show a dependency on the z coordinate, it is sufficient to plot the magnetic field in two dimensions, i.e., in the x-y plane with z = 0.

Positioning Performance-An Estimation-Theoretic Approach
Previously, we have calculated the magnetic field strength when placing elementary structures, in particular a cylinder and a sphere, in an outer homogeneous magnetic field.These magnetic fields H(a) depend on a variety of parameters, which we collect in a vector a = (a 1 , . . . ,a N ).Since the magnetic field is location dependent, the parameter vector a also contains location coordinates x.Thus, we can infer position information from the magnetic field strength measurement values.Usually, such measurements are noisy, which causes noisy position estimates.We are interested in the achievable performance of positioning based on noisy magnetic field strength measurements, in particular the variance in the noisy position estimates.
For the evaluation of the position estimation error, we apply the Cramér-Rao lower bound [23].The Cramér-Rao lower bound is a fundamental lower bound on the variance of any unbiased estimator.Its calculation in general requires the Fisher information matrix F with components where p Ĥ | a denotes the likelihood function with respect to the noisy magnetic field strength observations Ĥ, and E{•} defines the expectation operator with respect to that likelihood function.Finally, the variance of an estimate âi with respect to the parameter vector component a i is lower bounded by the ith diagonal element of the inverse Fisher information matrix.

Measurement Model
As already mentioned above, the magnetic field H depends on the position x.We are going to exploit this dependency for positioning.For that, we consider noisy measurements Ĥ = H(x) + ϵ (28) of the magnetic field strength, in particular its components.We assume that these measurements consist of the magnetic field strength H(x) itself, which is disturbed by additive white Gaussian noise ϵ = ϵ x , ϵ y , ϵ z T with zero mean and covariance In this case, the likelihood function is a conditional Gaussian probability density function, where n is the length of the error vector ϵ.

Fisher Information and the Cramér-Rao Lower Bound
With the spatial dependency of the magnetic field and the additive white Gaussian noise measurement model we have introduced in (28), the Fisher information matrix is calculated as [23] In the case of independent and identically distributed Gaussian noise, meaning Σ = σ 2 I, (31) simplifies to The Jacobian matrix contains the gradients of the magnetic field strength components with respect to the spatial parameters, which we wish to estimate.By applying the gradient operator for spherical coordinates to (17), we obtain the Jacobian matrix for a sphere in a homogeneous magnetic field, as introduced in Section 3.1.
For the cylinder case, we have calculated the magnetic field strength in Section 3.2.We observe that there is neither a dependency of the magnetic field strength on the height coordinate z nor a z component of the magnetic field strength itself.For this reason, position estimation in the z direction is impossible.So if we neglect the z dimension, we consider the remaining two-dimensional positioning problem in the x-y plane.By applying the gradient operator for the remaining polar coordinates (r, φ) to (25), we obtain the Jacobian matrix: With (32), the resulting Fisher information matrices are for the sphere and for the cylinder case with an identity matrix I 2×2 of dimensions 2 × 2. Inverting these Fisher information matrices results in the Cramér-Rao lower bound matrices for the sphere and for the cylinder.

Discussion of Results
Subsequently, we discuss the Cramér-Rao lower bound results of (38) and (39).The diagonal elements of the Cramér-Rao lower bound matrices provide the lower bounds for the variances of the parameter estimates.In our case, these parameter estimates are the position of the magnetic sensor in the particular coordinate system we have used.The corresponding variances describe the uncertainty along the unit vectors of the coordinate system.For the spherical coordinates, which we have used in (38), we obtain the variances in the radial (e r ), polar (e θ ), and azimuthal (e φ ) direction.In the case of the cylinder coordinates, used in (39), we obtain the variances in the radial (e r ) and azimuthal (e φ ) direction.Note, we omit the z coordinate for the cylinder case, since there is no dependency of the magnetic field on that coordinate.This makes an estimation of the z coordinate impossible, as already discussed in Section 4.2.
From both (38) and (39), we observe that the solutions factorize into parts that depend on the location of the sensor relative to the sphere with respect to the cylinder and parts that are location independent.We have formatted these factors to be unitless and therefore describe how the position estimation variances vary with respect to the squared radius of the sphere (r 2 sph ) or the cylinder (r 2 cyl ).Subsequently, we discuss these factors mainly in standard deviation form, i.e., the square root of these variance factors.

Dependency on Magnetic Sensor Properties
For both the sphere and the cylinder case, the solutions depend on the factor H 2 ∞ σ 2 .This factor does not depend on the sensor location and can be interpreted as a signal-to-noise power ratio for the magnetic field strength measurements, which we modeled in (28).This ratio is unitless and determined by the outer magnetic field strength and the measurement noise variance σ 2 , which is a measure for the sensor quality.As this signal-to-noise ratio increases by 10 dB, the position estimation variances decrease by one decade.In other words, the corresponding position estimation standard deviations decrease by 1 decade per 20 dB signal-to-noise ratio.

Dependency on Material Properties
A second location-independent factor contains the permeability µ r of the sphere's respective and cylinder's respective material relative to the outer medium.Figure 5 shows the material factor σ mat for both the sphere and the cylinder.In both cases, this factor shows a pole at µ r = 1.At µ r = 1, the permeability of the sphere's respective and cylinder's respective material is identical to that of the outer medium.Consequently, there is no distortion of the outer stimulating magnetic field.It remains homogeneous and therefore location independent, which makes position estimation impossible.This means that the estimation uncertainty in the form of the estimation variances is approaching infinity.For ferromagnetic materials, µ r ≫ 1, the material factor σ mat approaches an infimum of 1.As we can observe from (38) and (39), a relative permeability of µ r = 100 results in σ mat = 102 99 ≈ 1.03 for the sphere and σ mat = 101 99 ≈ 1.02 for the cylinder.So, there is only little to gain if the relative permeability is further increased.It is interesting to note that for ideal diamagnetic materials (µ r = 0), like superconductors for example, we observe different material factor values for the sphere (σ mat = 2) and the cylinder (σ mat = 1).

Dependency on Sensor Location
Previously, we have discussed those parts in the Cramér-Rao lower bounds that do not depend on the magnetic sensor location.Distortions in the magnetic field, caused by the sphere and the cylinder, are location dependent.From the Cramér-Rao lower bound theory, it is known that higher variations in the magnetic field strength provide better position estimation performance.Significant magnetic field strength variations can be observed around the magnetic structures.With increasing distance to those magnetic structures, the distortions relax, and the magnetic field becomes more and more homogeneous.Therefore, we expect higher position estimation performance around the magnetic object.Subsequently, we discuss the corresponding location-dependent factors in the Cramér-Rao lower bounds of (38) and (39).

Sphere
For the environment shown in Figure 1 with an outer homogeneous magnetic field in the z direction, we expect the resulting magnetic field to be rotationally symmetric with respect to the azimuth coordinate φ.From (13), we observe that there is neither a magnetic field component in the e φ direction nor a dependency on the azimuth coordinate φ of the remaining magnetic field components.Therefore, the position estimation performance is independent of the azimuth coordinate φ.
The radial dependency factor σ rad shown in (38) increases with the 4th power of the radial distance r normalized to r sph , which is the radius of the sphere.So, the position estimation error quickly increases when increasing the distance to the sphere.Figure 6 shows the radial dependency σ rad as a function of the normalized radial distance.The radial dependency factor affects the position estimation in all directions similarly.The diagonal elements of the Cramér-Rao lower bound matrix provide the lower bounds for the variances of position estimation errors in the radial, polar, and azimuthal directions, respectively.From (38), we observe that the corresponding standard deviation factors σ r , σ θ , and σ φ depend differently on the polar angle θ. Figure 7 shows the graphs of these factors.Additionally, we have plotted the magnitude of the position estimation error factors depending on the polar angle θ.We have normalized the polar angle-dependent diagonal elements in (38) to the magnitude's minimum.Therefore, with minima (min σ mag = 1) occurring at θ = 0 • and θ = 180 • .The performance factors are limited for the radial and polar directions.In contrast, the position estimation error in the azimuthal direction shows a singularity and goes to infinity for θ → 90 • .It is interesting to note that with σ φ = σ θ for θ = 0 • , 180 • .According to Figure 7, position estimation in the radial direction always performs better than in the polar or azimuthal directions.The block diagonal structure of the Cramér-Rao lower bound matrix indicates that the position estimation performance in the azimuthal direction does not depend on the position estimation performances achieved in the radial and polar directions.Figure 8 shows a three-dimensional sliced color plot for the magnitude of the locationdependent position error factor, which is defined as Magnetic flux lines indicate the magnetic field strength and orientation.In the slice plane y = 0, we observe the singularity discussed above.For the whole plane z = 0, or equivalently θ = 90 • , the factor σ φ and consequently σ loc go to infinity.Nevertheless, we still are able to estimate the remaining position coordinates, i.e., the radial distance r and the polar angle θ, with sufficient accuracy.

Cylinder
For the cylinder case, shown in Figure 3, the outer static magnetic field is directed along the x axis.The cylinder itself has infinite dilatation in the z direction.With this choice of the coordinate system, we obtain invariance with respect to the z coordinate and reduce the problem to a two-dimensional one.Consequently, position estimation in the z direction is not possible.
The radial dependency factor σ rad shown in (39) increases with the 3th power of the radial distance r normalized to the cylinder radius r cyl .Compared to the sphere, which has a finite size in all the three dimensions, the impact of the inhomogeneity on the outer magnetic field, and therefore the achievable position estimation accuracy, is more significant.Figure 6 shows the corresponding graph in comparison to the sphere case.
From (21) and Figure 4, we observe that the magnetic field shows both a magnetic field component in the e φ direction and a dependency on the azimuth coordinate φ of the magnetic field components.As indicated by the identity matrix in (39) however, the position estimation error variances are equal for both the radial (e r ) and the azimuthal (e φ ) directions.Despite the magnetic field not being rotationally symmetric, it is interesting to see that the position estimation performance is independent of the azimuth coordinate φ in this case.

Examples
Subsequently, we consider a vertical cylindrical pillar made of ferromagnetic material as an example and assess the achievable positioning accuracy when measuring the Earth's magnetic field in the area around that pillar.

Cramér-Rao Lower Bound in the Presence of Quantization Noise
We start from (39).The diagonal elements of the Cramér-Rao lower bound matrix are equal and provide the lower bounds for the variances of the position estimation in the radial (e r ) and azimuthal (e φ ) directions.We consider the square root of the diagonal elements and obtain as the Cramér-Rao lower bound for the estimation error standard deviation in both radial and azimuthal direction normalized to the radius of the cylindrical pillar.For the position estimation standard deviation of the Cramér-Rao lower bound, we obtain In our example, the outer magnetic field H ∞ is the Earth's magnetic field, which we consider as sufficiently homogeneous around the pillar.We may describe the signal-tonoise ratio (SNR) terms As discussed in Section 3.2, the ferromagnetic pillar only influences the magnetic field components perpendicular to its axis.So in our case, Next, we have a look at the measurement noise variance σ 2 H , or equivalently σ 2 B , which appears in the SNR terms in (45).As the source of measurement noise, we consider the quantization noise of the magnetic field sensor.The quantization noise variance can be calculated from the magnetic field sensor resolution as assuming that the quantization error is uniformly distributed within an interval of − B res 2 , + B res 2 .The range and resolution of several magnetic field sensors, which are currently used in smartphones, are reported in data sheets [24][25][26][27][28].The parameters relevant to us, together with the corresponding quantization noise variances, are summarized in Table 1.For completeness, we have included the SNR term B in Table 1. Figure 9 shows the Cramér-Rao lower bounds according to (44).The permeability of the cylinder's material is assumed to be sufficiently large such that the material factor becomes σ mat ≈ 1 with good approximation.
As an example, let us consider the sensor MMC246xMT, which provides a resolution of 25 nT.This resolution results in a quantization noise variance of 52.1 nT 2 and provides

dB with an outer horizontal
Earth magnetic flux density of B h = B ∞ = 21 µT.Note we express this variance in the physical unit 'Nanotesla'.We consider the physical unit 'Nanotesla' (nT) as a whole.To avoid misinterpretations in the notation, please note that a variance of 1 nT 2 = 1 (nT) 2 = 10 −18 (T) 2 .
To obtain a numerical example, let us assume a radius of r cyl = 10 cm for the vertical pillar.From Figure 9, we observe that the Cramér-Rao lower bound for the position estimation standard deviation at a distance of r = 10 • r cyl = 1 m is σ pos = 0.244 • r cyl = 2.44 cm.If we further increase the distance by a factor of ten, i.e., r = 10 m to the center of the cylinder, the position estimation standard deviation increases to σ pos = 24.4 m.
The position estimation standard deviation is proportional to the 3rd power of the (normalized) radius, thus meaning that the slope of the graphs in Figure 9 is three standard deviation decades per distance decade.For our example, position estimation at distances higher than 10 times the radius of the cylinder results in a standard deviation error, which rapidly exceeds the distance itself.
Note in Section 4 that we have derived the Cramér-Rao lower bounds for additive Gaussian noise, i.e., a Gaussian likelihood function.However, we have considered quantization and assumed a uniform likelihood function.In [29], it is shown that the Gaussian likelihood function provides the largest Cramér-Rao lower bound among all likelihood functions with equal variance.So, the Gaussian noise assumption provides a worst case.

Cramér-Rao Lower Bound in the Presence of Earth's Magnetic Field Noise
We continue our example of a vertical ferromagnetic pillar.Previously, we have assumed that the Earth's magnetic field strength, in particular its horizontal component, is known to us.The Earth's magnetic field, however, is fluctuating.We consider these fluctuations as stochastic processes.

Characterization of the Earth's Magnetic Field Noise
These stochastic processes are the sources of measurement noise in further investigations of our example.Figure 10 shows the fluctuations i.e., the differences with respect to the corresponding mean values, for the three components B x , B y , and B z of the Earth's magnetic flux density in a Cartesian coordinate system.The horizontal components B x and B y point toward the north and east, respectively.The vertical component B z points toward the Earth's center.The data was obtained from an observatory in Fuerstenfeldbruck, Germany, which is part of the INTERMAGNET network [30].The data was collected in May 2016, with a sampling interval of 1 min.Out of that data, we obtain the mean vector and covariance matrix of the horizontal components of the Earth's magnetic flux density for our further investigations.The variances of the fluctuations in the x and y directions are significantly different and show only little correlation.Another statistical property of the Earth's magnetic flux density that we are interested in is the shape of its probability density function (PDF).We compare the PDF of the Earth's magnetic flux density fluctuations with a corresponding Gaussian PDF with equal variance.For this comparison, we draw quantile-quantile (Q-Q) plots [31], as shown in Figure 11.If the sample data are Gaussian distributed, the Q-Q plot varies randomly around the main diagonal, which we have plotted as a dashed line.In the interval [−20, 20] nT, the quantiles of the samples are quite close to those of the corresponding Gaussian distribution, thus indicating a kind of 'bell' shape.Deviations indicate increased tails (x component) or skewness (z component) for instance.
We would like to emphasize that the analysis of the mean and variance of the Earth's magnetic flux density is based on an exemplary one-month data sample.Depending on the geomagnetic activity and geographical location, higher fluctuations may occur than in this example.Such higher fluctuations mean higher variance values, which are regarded as noise variance in our Cramér-Rao lower bound analysis.
To summarize, the PDFs of the Earth's magnetic flux density fluctuations are approximately, but not perfectly, Gaussian-shaped.However, as mentioned earlier, the Gaussian assumption provides a worst case assumption for calculating the Fisher information or, correspondingly, the Cramér-Rao lower bound.

Cramér-Rao Lower Bound
For the derivation of the magnetic strength field in Section 3.2, we have assumed an outer magnetic field in the x direction.The Earth's outer magnetic field shows components in both the x and y directions.With an appropriate rotation of the coordinate system, we can generalize the result in (25) and obtain with and σ mat , as defined in (39) and ( 48).Here, we describe the magnetic field in the form of the magnetic flux density B out = µ 0 H out instead of the magnetic strength H out .Both fields are proportional with the proportionality constant µ 0 .As a next step, we analyze how fluctuations of the magnetic field B ∞ transfer into magnetic field fluctuations in the presence of a ferromagnetic pillar.According to (49), we set B ∞ = B∞ + ∆B ∞ , insert this into (52), and obtain as noisy magnetic flux density measurements.We identify as the mean magnetic flux density and as the noise part.With (56) and (51), we can calculate the noise covariance matrix as For notational convenience, we have omitted the argument of matrix Φ.The Jacobian matrix of the mean magnetic flux density Bout is composed of two column vectors corresponding to the derivatives of Bout with respect to the radial and azimuthal directions.The derivative of matrix Φ is calculated as with ( 57) and (58), we can calculate the Fisher information matrix F = J T Σ −1 J according to (31).We obtain the Cramér-Rao lower bound matrix as with ∥ B∞ ∥ 2 = B2 x + B2 y as the squared magnitude of the mean magnetic flux density B∞ .The identities for the 2 × 2 matrices J −1 and J T −1 are derived in Appendix A.

Discussion of Results
The square root of the main diagonal elements of matrix C in (60) provide the Cramér-Rao lower bound for the standard deviation for the radial estimation error and for the azimuthal estimation error.For a more general description, the Cramér-Rao lower bounds in (61) and (62) are normalized to the radius r cyl of the cylinder.In Section 4, we have assumed a noise covariance of form Σ = σ 2 I.With this assumption, the noise variance is constant, i.e., independent of the location of the sensor, and its components in both the x and y directions are equal.However, the noise covariance, shown in (57) for our example, depends on the location of the sensor through variable r and the argument of matrix Φ, which in general show different main diagonal elements through (51).Compared to the evaluations in Section 4, the Cramér-Rao lower bounds for both the radial and azimuthal estimations are dependent on the azimuth φ as well.Figure 12 shows the Cramér-Rao lower bounds for the position estimation standard deviations of the radial and azimuthal components for a vertical cylindrical pillar in the Earth's magnetic field according to (61) and (62) for our example.The dependency on the azimuth φ is clearly visible, since the contour lines are not circularly shaped.For the calculation of the position estimation standard deviation of the Cramér-Rao lower bound we have started with the 2 nd line in (60) and used the identity tr J −1 ΣJ = tr ΣJJ −1 = tr(Σ) for the trace operator.As in the case of the radial and azimuthal estimations, the position estimation of the Cramér-Rao lower bound according to (63) depends on the azimuth φ through matrix Φ.However, the influence of matrix Φ for the calculation of the Earth's magnetic field noise covariance matrix Σ rapidly decreases with the 2nd power of radius r, i.e., the distance of the sensor to the center of the cylinder.Figure 13 shows the position estimation of the Cramér-Rao lower bound according to (63).A significant dependency on the azimuth φ is not visible.For sufficiently large distances r, the position estimation of the Cramér-Rao lower bound approaches the asymptotic value which is similar to (45) in its structure.In particular, we observe that the asymptotic position estimation of the Cramér-Rao lower bound is independent of the azimuth φ.For our example, the SNR term in (64) becomes Note that the fluctuations in the Earth's magnetic flux density are accounted for by the covariance matrix Σ ∞ .Higher fluctuations lead to a higher trace, tr(Σ ∞ ), of that covariance matrix and consequently to a lower SNR and an increased Cramér-Rao lower bound.

Cramér-Rao Lower Bound in the Presence of Quantization and Magnetic Field Noise
For the consideration of both the quantization noise and the Earth's magnetic field noise, we calculate a general noise covariance which we then can insert into (61)-( 63), thus replacing Σ ∞ .With the addition of the corresponding covariance matrices in (66), we consider the magnetic field noise and sensor quantization noise to be uncorrelated.For the asymptotic case, we have plotted the Cramér-Rao lower bound according to (64) but replaced Σ ∞ with Σ in (66).The graphs are shown in Figure 14 as solid lines.For comparison, we have included the results that have been obtained for the different sensors in Section 5.1 (Figure 9), where only the quantization noise has been taken into account (dashed line with markers).When comparing the magnetic field noise variance tr(Σ ∞ ) = 570.36nT 2 with the sensor quantization noise tr σ 2 B summarized in Table 1, we observe that for the sensors AKM8975, YAS529, and YAS532, the sensor quantization noise is dominant.Therefore, we approximately obtain the performance derived in Section 5.1 for the sensors AKM8975, YAS529, and YAS532.This is clearly visible in Figure 14.The results when considering quantization and magnetic field noise match with the corresponding results when taking into account quantization noise only (black, magenta, and blue dashed lines with markers), and the graphs overlap.For the sensor HMC5983, the Earth's magnetic field noise and the sensor quantization noise variances are within a similar order of magnitude.An additional performance degradation when taking into account the Earth's magnetic field noise (red line) in addition to the sensor quantization noise can be observed in Figure 14.When applying sensors with significantly lower sensor noise, like the MMC246xMT, the Earth's magnetic field noise becomes dominant.Therefore, we can expect the results derived in Section 5.2 to have good approximation.Cramér-Rao lower bound for position estimation standard deviation for a vertical cylindrical pillar in the Earth's magnetic field according to (64), with (66) and σ mat = 1.Solid lines represent results for considering the Earth's magnetic field noise and the sensors' quantization noise.Dashed lines with markers ('+') repeat results from Figure 9, where only quantization noise has been considered.In our evaluations above, we have considered sensor quantization and fluctuations in the magnetic flux density of the Earth as sources of noise.Electrical currents in the environment can also lead to impairments in magnetic field-based position determination.The resulting magnetic fields are superimposed on the Earth's magnetic field and lead to higher magnetic field noise and therefore to higher CRLBs.

Summary and Conclusions
In this paper, we aimed at assessing the achievable positioning performance based on measuring the ambient magnetic field around elementary ferromagnetic structures, in particular a sphere and a cylinder, from a theoretical point of view.For that, we have calculated the Cramér-Rao lower bound for the position estimation around a sphere and a cylinder placed in an outer homogeneous magnetic field.For these structures, the magnetic field is analytically known.
We have seen that usual ferromagnetic materials like ordinary steel with a relative permeability of several thousands are sufficient.Using specifically designed ferromagnetic materials like permalloy, a nickel-iron magnetic alloy with µ r ≈ 10 5 , provides only little additional performance improvements.
The results have also shown that the position estimation error variance is proportional to the measurement noise variance.It is obvious that the sensor itself is a source of measurement noise.In this paper, we have considered the sensor quantization noise, which can differ by several orders of magnitude for different sensors.For the case that the instantaneous outer magnetic field strength is not known perfectly, we have considered outer magnetic field fluctuations as another source of measurement noise.In an example, we have observed an outer Earth magnetic field variance of 570.36 nT 2 of the horizontal component.For low resolution sensors, the outer magnetic field variance can be neglected.For high resolution sensors, the outer Earth magnetic field might be dominating.In order to take advantage of a high sensor resolution, the instantaneous Earth magnetic field strength has to be known.A potential solution is to simultaneously estimate that value with the help of a magnetic sensor array instead of a single sensor.
The influence of a magnetic structure to the ambient magnetic field relaxes with increasing distance.Therefore, the position estimation performance also decreases.For a ferromagnetic sphere in an outer homogeneous magnetic field, we have seen that the position estimation error variance increases by eight decades per distance decade.Correspondingly, the position estimation standard deviation increases by four decades per distance decade.Compared to a sphere, which has a limited size in all three dimensions, we additionally investigated a cylinder with infinite size in one dimension.For the cylinder case, the position estimation standard deviation increases by three decades per distance decade.However, position estimation is only possible in two dimensions, i.e., the plane perpendicular to the cylinder.Due to the rapidly increasing error, position estimation with sufficient accuracy is possible only in close proximity of such ferromagnetic structures.Both the position estimation standard deviation and the corresponding distance are proportional to the size of the structure, which is in our case its radius.
In this paper, we have focused on elementary structures that allow for the analytic description of the ambient magnetic field.This provided us with some general insights into the positioning error and its principle dependencies on relevant parameters like the permeability of the material.In practice, however, we face structures that are much more complex.Such structures require a numerical calculation of their ambient magnetic field and are out of the scope of this paper.Once the magnetic field has been calculated, the methods used in this paper for assessing the achievable position estimation performance, in particular the Cramér-Rao lower bound, can be applied in a similar way.Undoubtedly, it is interesting to experimentally verify the research questions raised in Section 1.2 for the considered setups of a ferromagnetic sphere or cylinder in an outer homogeneous magnetic field.As we focus on the theoretical aspects in this paper, such experimental support is the subject of further work.

Figure 1 .
Figure 1.Solid sphere with permeability µ in in an outer medium with permeability µ out .

Figure 2 .
Figure 2. Magnetic field strength and magnetic flux lines for a sphere with relative permeability µ r = 100 in an outer homogeneous magnetic field.

Figure 3 .
Figure 3. Solid cylinder with permeability µ in in an outer medium with permeability µ out .

Figure 4 .
Figure 4. Magnetic field strength and magnetic flux lines for a cylinder with relative permeability µ r = 100 in an outer homogeneous magnetic field.

Figure 5 .
Figure 5. Dependency of position estimation standard deviation on the relative permeability µ r .

Figure 6 .
Figure 6.Radial dependency of position estimation standard deviation.

Figure 7 .
Figure 7. Polar angle dependency of position estimation standard deviation for a sphere in a homogeneous magnetic field.

Figure 8 .
Figure 8. Location dependency σ loc of position estimation standard deviation for a sphere in a homogeneous magnetic field.

H 2 ∞/2σ 2 H = B 2 ∞/2σ 2 B
in (45) either in the form of the magnetic strength H or the magnetic flux density B. Both fields are proportional, where the proportionality constant is µ 0 = 4π • 10 −7 V s A m in a vacuum, which equals the permeability of air with good approximation.The Earth's magnetic flux density in Central Europe is approximately B = 48.3µT, with a horizontal component of B h ≈ 21 µT and a vertical component of B z ≈ 43.5 µT.

Figure 9 .
Figure 9. Cramér-Rao lower bound for position estimation standard deviation versus the distance of the sensor to the center of a vertical cylindrical pillar in the Earth's magnetic field according to (45) with σ mat = 1.
B [nT] B x − Bx B y − By B z − Bz

Figure 11 .
Figure 11.Quantile-quantile plots for the Earth magnetic field fluctuation in Fuerstenfeldbruck, Germany during May 2016.

Figure 12 .
Figure 12.Cramér-Rao lower bound for position estimation standard deviation of the radial and azimuthal component for a vertical cylindrical pillar in the Earth's magnetic field according to (61) and (62), with σ mat = 1.(a) Radial estimation; (b) Azimuthal estimation.

Figure 13 .
Figure 13.Cramér-Rao lower bound for position estimation standard deviation for a vertical cylindrical pillar in the Earth's magnetic field according to (63) with σ mat = 1.
Figure14.Cramér-Rao lower bound for position estimation standard deviation for a vertical cylindrical pillar in the Earth's magnetic field according to (64), with (66) and σ mat = 1.Solid lines represent results for considering the Earth's magnetic field noise and the sensors' quantization noise.Dashed lines with markers ('+') repeat results from Figure9, where only quantization noise has been considered.