Summary
Magnetic field localization is based on the fact that the Earth’s magnetic field is distorted in the vicinity of ferromagnetic objects. When ferromagnetic objects are in fixed positions, the distortions are also fixed and, thus, contain location information. In our prior work, we proposed a simultaneous localization and calibration (SLAC) algorithm based on a Rao-Blackwellized particle filter that enables magnetic train localization using only uncalibrated magnetometer measurements. In this paper, a lower-complexity version of the SLAC algorithm is proposed that only estimates a subset of calibration parameters. An evaluation compares the full and reduced SLAC approach to a particle filter in which the magnetometer is pre-calibrated with a fixed set of parameters. The results show a clear advantage for both SLAC approaches and that the SLAC algorithm with a reduced set of calibration parameters achieves the same performance as the one with a full set of parameters.
1 INTRODUCTION
In current railway systems, the safety distance between two trains is at least the absolute braking distance at all times. High-speed trains have to be separated by several kilometers due to the high mass of the vehicles. These large headways lead to inefficient use of the track network and reduce the overall capacity. While this may be tolerable in rural areas, it is not in urban areas in which the demand for transport capacity is high. To address this issue, network capacity has to be increased either by using the existing tracks more efficiently or by building new tracks. The latter is often not feasible because of the costs and required space, which is scarce in urban areas. It is, therefore, desirable to increase capacity by reducing the distance between trains below the absolute braking distance. The reduction of safety distances requires accurate real-time position information from all trains in a certain area to maintain traffic safety. In contrast to road vehicles, where it is usually sufficient to observe the environment of the vehicle, (e.g., with cameras or radars) to maintain safe operation, trains rely on information beyond the measurement range of such sensors due to large braking distances. Furthermore, direct distance measurements between trains may result in distances shorter than the actual distance due to the track geometry. Nowadays, train localization is based on the deployment of costly dedicated infrastructure. It is, therefore, desired to have a low-cost infrastructure-less alternative. One more modern alternative is the use of global navigation satellite systems (GNSSs) like GPS, Galileo, and GLONASS (Marais et al., 2017). While GNSS is a good choice in many environments, it is also vulnerable to jamming, spoofing, and shadowing. Fortunately, jamming and spoofing are manmade and encountered only rarely. In contrast, shadowing is encountered regularly and leads to reduced performance and availability. In the extreme case of shadowing (e.g., in tunnels), GNSS signals are completely blocked. But even in urban canyons, the amount of visible satellites and the quality of their geometry are already dramatically reduced. When shadowing is encountered only for short periods, GNSS can be aided with an inertial sensor or an odometer. Depending on the quality of the sensors, it is possible to maintain a certain position accuracy during shadowing for a couple of seconds up to a few minutes. For longer outages, costly sensors with a very high quality are required.
In our works, Heirich et al. (2017) and Siebler et al. (2020), we proposed to exploit position-dependent distortions of the Earth’s magnetic field for train localization to complement GNSS. These magnetic distortions are a result of ferromagnetic material in the existing railway infrastructure and are free to use. The use of such distortions was also proposed for indoor (Frassl et al., 2013; Haverinen & Kemppainen, 2009; Jung et al., 2015; Kok & Solin, 2018; Li et al., 2012; Shao et al., 2016), car (Shockley & Raquet, 2014), and airplane (Canciani & Raquet, 2017) localization.
Assuming a map of the magnetic distortions exists for the whole track network, it is possible to localize a train on a track by matching the measurements of a low-cost onboard magnetometer with the magnetic map. One of the issues encountered during our research is that the matching of the magnetometer data to a map requires a calibration of the sensors. This is because the magnetometer is influenced by ferromagnetic material in the vicinity of its mounting position. A lack of calibration causes the measurements of sensors in different trains to differ even when they measure along the same trajectory. Sensor calibration is, therefore, crucial if the same map should be used for different trains. Common magnetometer calibration procedures require the rotation of the sensor and the platform it is mounted on in a homogenous magnetic field (Kok & Schön, 2016; Renaudin et al., 2010; Vasconcelos et al., 2011; Wahdan et al., 2015). The rotation ensures that all sensor axes are excited and the calibration parameters become observable. Alternatively to rotating the sensor, the sensor can be exposed to a known but varying magnetic field. Due to the high mass and the large volume of the train, both approaches require considerable effort. In Siebler et al. (2021a, 2021b), we showed that the magnetic distortions along a railway track can be used to render the calibration parameters observable. This requires a map of the magnetic field along the track where the map itself can be recorded with an uncalibrated sensor. During calibration, the current train position is used to obtain the magnetic field from the map that the magnetometer on the train should measure. Since we are mainly interested in estimating the train position, we cannot assume the position to be known. Therefore, simultaneous localization and calibration (SLAC) is performed by a Rao-Blackwellized particle filter. In this paper, a reduced version of the SLAC algorithm that only estimates a subset of the calibration parameters is introduced. In an evaluation, the different variants of the algorithm are compared to a particle filter that uses a fixed set of calibration parameters.
2 METHODOLOGY
2.1 Frame Definitions
For the subsequent presented methodology, we first define the used coordinate frames. All frames are attached to the sensor but have different orientations.
Track frame: This frame is always parallel to the track plane. The x-axis points into the direction of an increasing along-track position, s, the y-axis to the right, and the z-axis points downwards.
Vehicle frame: The x- and y-axes are parallel to the floor of the train. The x-axis points along the train, the y-axis to the right, and the z-axis points downwards. Due to the definition of the frames, for zero roll and pitch angles between the track and vehicle frame, the x-axis of the vehicle frame is either pointing in the same or opposite direction as the x-axis of the track frame. The same holds for the y-axis. Therefore, the yaw angle is limited to 0° or 180°. The actual yaw angle depends on how the train is placed on the track, which is information that is typically known or easy to obtain. To simplify notation, in the remainder of the paper, the vehicle frame is always assumed to be rotated such that the yaw angle to the track frame is zero.
Body frame: This frame coincides with the physical x-, y-, and z-axes of the magnetometer triad. The relation between the body and vehicle frame is assumed to be fixed and defined by the mounting attitude of the sensor in the train.
The frame definitions are also visualized in Figure 1. Regarding the notation, a superscript on a vector indicates in which frame it is measured (e.g., the variables at and av define the same measurement a in the track frame t and the vehicle frame v). A vector is transformed from one frame to another by the rotation matrix between those two frames. For the previous example, we can write 
and 
.
(left) track frame; (right) vehicle frame
2.2 Localization Particle Filter
Magnetic field-based localization utilizes time persistent spatial variations in the Earth’s magnetic field. In its essence, the localization can be seen as a pattern-matching problem. The goal is to find a sequence of positions on the track for which the measurements of a train-mounted magnetometer best fit the magnetic field stored in a prior recorded map. In this paper, the localization problem is addressed by a Bayesian filter that attempts to calculate the posterior probability density function (pdf), p(d0:k | z1:k, m), of a state vector, d, from time step 0 to k given the magnetometer measurements, z, of the magnetic flux density vector and a magnetic map, m(·). For simplicity, we assume for now that the measurements are measured in the same frame as the map and that the sensor is calibrated. Under the given assumptions, a normal particle filter, e.g., a sampling importance resampling filter (Arulampalam et al., 2002), is suitable to estimate the posterior density with a weighted particle set, 
:
1
where Np is the number of particles, 
the weight associated to particle i at time step k, and 
is the Dirac distribution, which vanishes everywhere in the state space except at the value of particle i. The filter algorithm is composed of two steps. First, in a prediction step, new samples are drawn from a proposal distribution. Here, the chosen proposal distribution is equal to the movement model of the train and hence:
2
Second, for each predicted particle, the weight is updated based on the newest measurement by multiplying the old weight with the measurement likelihood:
3
For a calibrated sensor, the measurement model depends only on the current state:
4
where zk is the magnetic flux density vector measured by a magnetometer triad and 
is a white Gaussian measurement noise. Thus, the likelihood depends only on the current state and the weight update reduces to 
. With Equation (4), the likelihood becomes 
. The likelihood depends on the Mahalanobis distance, 
, between the magnetic vector in the map at the particle position and the current measurement. Therefore, the likelihood is high for particles at a position where the map is close to the measured magnetic vector and vice-versa. This shows that sensor calibration is crucial for the filter to work. For example, in the presence of a bias or scale errors, the distance is not guaranteed to be small when the particle is in the correct position.
2.3 Magnetometer Calibration
The magnetometer under consideration in this paper has three orthogonal axes. This enables the sensor to measure the direction and amplitude of the magnetic field at the current position. One of the difficulties in measuring the magnetic field is caused by the disturbance from nearby magnetic material. Hard iron effects will cause an offset in the measured magnetic field vector and soft iron effects interact with the external magnetic field, resulting in a rotation and scaling of the measured values. These two effects can be captured with a linear transformation of the true magnetic field (Kok & Schön, 2016):
5
where 
is the magnetic flux density vector measured in the sensor’s body frame at the discrete time step k, 
is the true magnetic field in the track frame, and 
is the rotation matrix from the track to the body frame. The calibration parameters are contained in a matrix, 
, that accounts for soft iron effects and a vector, 
, accounting for hard iron effects. As in Equation (4), nk is the measurement noise of the magnetometer. Transforming the map from the track to the body frame with a single rotation implicitly introduces certain assumptions. First, for Equation (5) to hold, the sensor during the map creation has to be mounted at the same height and lateral position as the sensor that is measuring 
. Second, the sensor position is not affected by the rotation. The first assumption can be taken care of by placing the sensors always at the same position relative to the track, for example, in the center of the train’s under-floor. The second assumption is never fulfilled completely. If the train has a non-zero roll angle, the sensor will move away from the track center to the left and to the right. Fortunately, the roll angle of a train is typically small, e.g., in the example used in the evaluation, the angle was below ≈ 2.5°. For a sensor mounted on the floor at a one-meter height, the change in the sensor position is only a few centimeters. In comparison, to observe a significant change in the magnetic field with regard to the measurement errors, the position has to change in the order of meters. Therefore, the change of the sensor position due to rotations is ignored.
The sensor model in Equation (5) is linear in its parameters and calibration can be performed with a linear least-squares estimator assuming that the true magnetic field, 
, and the rotation to the body frame are known. This becomes clearer when the model is rewritten as a function of a single parameter vector comprised of the elements of C and b:
6
where 
is the value a calibrated sensor should measure and H is a matrix:
7
and 
is the parameter vector:
8
where ci: is the i-th row of matrix C and bi is the i-th component of the bias vector. Equation (6) is clearly linear in θ, and therefore, can be estimated by a linear least-squares estimator. The least-squares estimator is obtained by stacking the left-hand and right-hand sides of Equation (6) for Nz measurements:
9
and then calculating the pseudo-inverse, (DT D)−1DT, of the design matrix, D. If the inverse (DT D)−1 exists, the least-squares solution is unique. For the example we investigate in this paper, this is the case from which we conclude that the magnetic variations contain enough information to fully observe the calibration parameters. Because the parameter vector has dimension 12 and a single measurement in the sensor’s body frame, 
, has dimension three, at least four different measurements are required. For a good estimate, in practice, we will use more measurements to mitigate the errors introduced by the measurement noise. Furthermore, for the least-squares estimator, the calibration parameters are considered constant. If this is not the case, the Kalman filter can be used to track the calibration over time. The implementation only requires a state space model that describes how the parameters change over time.
2.4 Simultaneous Localization and Calibration (SLAC)
In this section, the SLAC algorithm is developed and it is shown that, in the railway domain for calibration, it is not necessary to record the map with a calibrated magnetometer.
2.4.1 Mapping and Calibration
The main problem encountered during the calibration of a train-mounted magnetometer is obtaining the true magnetic field in the body frame, 
, a calibrated sensor is supposed to measure because the rotation matrix, 
, and the true magnetic field, 
, in the track frame can continuously change along a railway track. Hence, to obtain the correct value, three quantities are required:
Along-track position sk of the magnetometer on the track at time step k
Map function mt (sk) of the magnetic field along the railway track that returns the true magnetic field
from Equation (5) in the track frame at position skRotation matrix
between the body frame and the track frame of the sensor at time step k
The map creation is a problem for itself, particularly in GNSS-denied areas such as tunnels. One possibility in such areas is to equip a train with additional, possibly costly, sensors to perform simultaneous localization and mapping (SLAM). The map creation, therefore, will require some effort, but the calibration approach in this section makes it possible to use the same map for different train types. In addition to the position, to create the true magnetic field map in the track frame mt (sk), a calibrated magnetometer and its attitude, 
, relative to the track are required. That means the mapping of the true magnetic field requires a calibrated magnetometer and, to calibrate the magnetometer, the true magnetic map is required. To resolve this issue, we show in the remainder of this section that, for the railway scenario instead of the true map, mt (·), it is sufficient to record the map with an uncalibrated sensor and then perform SLAC with respect to this map. In order to prove this, in a first step, we define the uncalibrated magnetic map, 
, in terms of the true magnetic map, mt (·), and the calibration parameters, 
and 
, of the magnetometer used to create the map:
10
where 
is the position-dependent rotation matrix between the track frame and the body frame of the magnetometer during map creation. To differentiate between the frames of the magnetometer that is to be localized and the frames of the magnetometer during mapping, the latter are indicated by 
and 
. Note, Equation (10) is closely related to Equation (5), but instead of giving the relation between the measurement at time step k and the true magnetic field, Equation (10) shows the relation between the uncalibrated map in the body frame and the true map in the track frame at position sk.
In a second step, the link between the measurements used for localization, 
, and the uncalibrated map is created by solving Equation (10) for mt (sk) and plugging the result for 
in the sensor model in Equation (5), which gives the relation between the uncalibrated map, 
, and the measurements used for localization:
11
From Equation (11), one can see that the mapping between the map and the measurement is, again, linear but time variant. Fortunately, a train is basically aligned with the pitch and yaw angle of the track and the only major change in rotation that might happen on the roll angle depends on the train type and speed. From our experience, this difference is also rather small, and we will show in the evaluation that, for the evaluated data set, the difference is below 0.5°. Hence the v and 
frame coincide and, with 
, we can write:
12
which is now a linear model with constant parameters. Because the map 
in Equation (12) is in the body frame of the sensor during mapping, the map not only contains magnetic variations due to inhomogeneities, but also the variations due to a changing attitude. Luckily, the assumption 
enables the use of these attitude-dependent variations without having to estimate the platform attitude.
The simplification 
is particular for railways. For example, when the magnetometer is mounted in a road vehicle, the attitude of the vehicle frame with respect to the road during mapping and localization can change more freely and, hence, the rotations will not cancel each other out. In such cases, it is necessary to estimate the rotation matrix alongside the position and calibration parameters. Furthermore, the map should be recorded with a calibrated sensor and with a known attitude to create a map of the true magnetic field. Otherwise, the resulting calibration parameters will depend on the sensor’s attitude during mapping and localization. In our future research, we will investigate the feasibility of SLAC for other platforms than trains.
2.4.2 SLAC With Rao-Blackwellized Particle Filter
So far, the train position sk was assumed to be known. Since the primary goal of this paper is to estimate the train position based on measurements of an uncalibrated magnetometer, this is not realistic. This means the problem we are facing is that, for position estimation, the magnetometer has to be calibrated and, to calibrate the magnetometer, the train position is required. This results in a chicken-or-the-egg problem. In Siebler et al. (2021b), we presented an approach to this problem that will be briefly introduced. The approach is based on the observation that the SLAC problem is related to the SLAM problem. In SLAM, the goal is to estimate the static positions of landmarks based only on the observations of a robot that, at the same time, uses the landmark positions to localize itself. In SLAC, we now replace the unknown static landmarks’ positions with the calibration parameters. Therefore, we can employ the idea of FastSLAM (Montemerlo et al., 2002). In FastSLAM, the joint posterior of the robot and landmark positions is decomposed into two parts and, then, estimated with a Rao-Blackwellized particle filter. In the SLAC approach, we do the same and decompose the joint posterior of the dynamic train state d from Equation (1) and the calibration vector θ from Equation (8) to obtain:
13
The first term on the right side of Equation (13) is the posterior density of the parameter vector, θk, given the history of the train state, the magnetometer measurements up to the current time step k, and the map. In the posterior of the parameter vector, the train state containing the position is assumed to be known, and the parameter vector is allowed to be time variant as indicated by the time index k. According to Section 2.3, the posterior therefore is obtainable with a linear Kalman filter. The state space model for the Kalman filter assumes a random walk for the parameters:
14
with process noise 
and 
. The process noise covariance is a tuning parameter and influences how fast the parameters can change. The measurement model is in the form of Equation (6), where 
is replaced with the values in the map at the current train position:
15
The second term in Equation (13) is the posterior of the train state given the measurements and the map. Due to the nonlinear relation between the magnetic field and the train position here, a particle filter as described in Section 2.2 is utilized that approximates the posterior by a weighted set of particles:
16
In the particle filter update step in Equation (3), the old weight is multiplied with the likelihood of the measurement given the train position and the map:
17
where 
is the state trajectory of particle i. After some thought, one might notice that the likelihood in the weight update is not conditioned on the parameters but, without them, the likelihood given by the measurement model in Equation (15) is not defined. Hence, we need a way to incorporate the parameters into the likelihood without assuming them to be known. This can be achieved with marginalization over the estimated calibration parameters:
18
in the marginalization is the predicted Gaussian pdf of a Kalman filter conditioned on the trajectory 
of particle i:
19
with the estimated posterior mean, 
, and covariance, 
, at time step k–1. Note that both the estimated parameters and state covariance matrix depend on the particle trajectory due to the measurement matrix H(·). The integral in Equation (18) is again Gaussian because both pdfs are Gaussian and the relation between 
and θk is linear, thus, the likelihood becomes:
20
with 
. The marginalization adds the uncertainty of the calibration parameters to the uncertainty due to measurement noise. The likelihood, therefore, becomes narrower when the parameters are estimated with a high certainty and wider otherwise. The marginalization of the likelihood in the update step is done also in FastSLAM, but for SLAM, the marginalization is with regard to the estimated landmark positions rather than the calibration parameters.
From an implementation point of view for each particle, a separate 12-dimensional Kalman filter is required that tracks the corresponding mean and covariance of the parameters required in the particles’ weight update. To reduce computational complexity, we exploit that the state vector of the Kalman filter can be decomposed into three independent substates with dimension four due to the chosen state space model of the parameters. Each substate is associated with only one element of the measured magnetic vector, 
, and contains four parameters: one row of the calibration matrix 
, and one element of the bias vector 
. For each substate, a separate Kalman filter per particle is required. Due to the quadratic complexity of the Kalman filter with respect to the state dimension, three separate filters lead to a reduced overall complexity.
Algorithm 1 shows the high-level pseudocode of SLAC where the Kalman gain matrix, Kk, is given by 
and 
is the innovation covariance.
SLAC
2.4.3 Reduced SLAC
Further complexity reduction can be achieved by reducing the numbers of estimated parameters. When the soft-iron effects and misalignment errors are small compared to the scale and bias errors, one can only estimate a scale and bias for each axis during calibration. Hence, the estimated state vector becomes:
21
and the measurement matrix reduces to:
22
with the i-th element of the true magnetic field in the body frame, 
, and the diagonal elements of the calibration matrix, cii. When soft-iron and misalignment effects are negligible, this implies a diagonal structure of 
and C. In the SLAC approach, this is not yet helpful because the matrix 
is estimated, which is a product of the calibration matrices of two different magnetometers and their relative attitudes. Therefore, we need to have a closer look at the conditions under which 
also becomes diagonal. The simplest example is when the different sensors are mounted with the same orientation so that 
and 
cancel each other out. If this is not the case, one can try to rotate 
into the body frame of the map by left-multiplying Equation (12) with the inverse of 
to obtain the rotated measurement vector:
23
Because 
and C are assumed to be diagonal, 
is diagonal when the similarity transform 
preserves the diagonal structure of C. This is fulfilled only for some rotation matrices and, hence, it is recommended to mount the different magnetometers in the same orientation. Due to the small size of low-cost magnetometers, this should not be a problem.
Fortunately, in the example data set presented in the next section, the rotation between the two body frames was given by the rotation matrix:
24
For this rotation, the similarity transform preserves the diagonal structure of C:
25
Considering Equations (12) and (25) as well as a diagonal 
, it follows that also 
is diagonal, and only six parameters have to be estimated in the Kalman filter of the Rao-Blackwellized particle filter. As before, the state is split into three parts: one for each axis. This results in three Kalman filters with dimension two per particle. For reduced SLAC, the pseudocode in Algorithm 1 needs only two changes. The reduced measurement matrix in Equation (22) is used instead of Equation (7), and when the new measurement is obtained, it is immediately rotated into the map frame by left-multiplication with the matrix 
.
2.4.4 Applicability of the Approach
For SLAC to be applicable, it has to be ensured that the assumption 
is fulfilled and that the sensors are mounted roughly at the same height. This ensures that the different trains measure the magnetic field at the same localization. When it comes to the train type, we would expect that SLAC works for different engines types. For electrical engines, one might have to be careful when it comes to the sensor mounting positions. The sensor should be mounted at a position at which it is not affected by the power electronics of the engine. In addition, in our prior work, Siebler et al. (2020), we saw that magnetometers measure a changing magnetic field due to the current in the overhead lines. Fortunately, these currents have known frequencies that can be filtered out. It should be noted that some trains are powered over a DC power line. In such cases, the current induced magnetic fields are not so easy to filter out, rendering it problematic to perform SLAC. Besides different engine types, trains can also have largely different dimensions but it is expected that this has no big effect as the calibration parameters are affected only by materials in the vicinity of the sensor.
3 EVALUATION
In this section, the proposed algorithms are evaluated with measurements from two rural trains. The results of SLAC and reduced SLAC are also compared against a simple particle filter that ignores the calibration and a particle filter that uses a fixed set of calibration parameters.
3.1 Measurement Campaign
The data for the evaluation was collected in a measurement campaign performed on the track network of the Harzer Schmalspurbahnen. During the campaign, magnetic measurements were recorded with the trains shown in Figure 2. Both trains were equipped with a single-frequency u-blox LEA-M8T GNSS receiver and a magnetometer triad contained in an Xsens MTi-G700 inertial measurement unit (IMU). The magnetometer on the steam locomotive was mounted on the back in a steel cabinet. In the railcar, the magnetometer was mounted on the floor. The sensors were placed in the center of the track and had a height above the rails of roughly 1.5 m in the steam locomotive and 1 m in the railcar. The relative attitude between the body frames of the sensors during the steam locomotive and railcar measurements are given by Equation (24), which corresponds to a pitch angle of 90° (in the steam locomotive, the sensor is hanging on the back of the train and, in the railcar, the same sensor is placed flat on the floor). The railcar and the steam locomotive traveled at speeds ≤ 50 kmh−1. The evaluation is based on a 7.2 km long track segment between Quedlinburg and Gernrode. A satellite image of the track is shown in Figure 3. During the recording of the localization data, the railcar was accelerating to a speed of ≈ 40 kmh−1. The different filters were started the moment the train started driving. This is important for SLAC because, in standstill, the parameters are not observable and, due to the system model in the Kalman filter of Equation (14), their uncertainty keeps growing in standstill which could potentially lead to a divergence of the filter. The same happens when the train comes to a stop. The simplest approach to avoid this issue is to start and stop the filter in such situations.
(top) steam locomotive with magnetometer mounted on the backside inside a steel cabinet, as shown in the right part of the picture; (bottom) diesel railcar
The track in the evaluation is ≈7.2-km long and runs through a rural area between Quedlinburg and Gernrode with a mostly open-sky view (Map data: Google Earth, GeoBasis-DE/BKG).
3.2 Magnetic Field Map
The magnetic vector field in the map is based on the magnetometer data from the steam locomotive. During the map creation, the ground truth was obtained from the GNSS receiver. In the map creation, the GNSS speed was integrated to transform the magnetometer measurement time series into a series which each element was a 1D along-track position on the track assigned to it. This assignment transforms the measurement series from the time into the spatial domain. Based on the transformed time series, noise was removed by fitting a linear combination of Gaussian basis functions to it. The basis functions are placed equidistant along the track with a fixed bandwidth. This procedure ensures a smooth map and gets rid of small-scale spatial fluctuations in the measurement series.
Quick access to the map is critical for the overall performance because, in every update step, the magnetic field at all particle positions is required. Instead of accessing the continuous map function defined by the basis functions, the map is discretized at an equidistant grid with a spacing of Δs = 10 cm. The discretization results in an array that contains the magnetic vector for each (discrete) along-track position on the track. Access to the map is performed by indexing. The index, Ik, for a certain position, sk, is easily calculated by Ik = round(sk / Δs) +1, assuming the first array index is one. For the evaluation, the map additionally contains the corresponding Earth-centered, Earth-fixed (ECEF) position from the GNSS receiver.
3.3 State Space Model and Parametrization
The train state is modeled as a discrete Wiener process acceleration movement model:
26
where 
is white Gaussian noise with covariance matrix:
27
The parameter σw is set to approximately the maximum acceleration change within sampling period T, as suggested in Bar-Shalom et al. (2002). For a train, the jerk 
(change of acceleration) is usually limited to a value of 1 ms−3. The maximum change in the acceleration during one sampling period, therefore, is σw = T · 1 ms−3. For the system model of the parameters in Equation (14), the covariance matrix is a tuning parameter and was set experimentally to σu = 0.005, where the dimension of the identity matrix I depends on the number of estimated parameters.
The particles for all filters are randomly initialized. The samples for the different particles and state components are drawn independently from uniform distributions centered at the corresponding true state. For the position, the width of the uniform distribution is 20 m for a speed of 1 ms−1 and an acceleration of 0.5 ms−2. The Kalman filter initial state covariance for the parameters is set to one for all states. In the initial Kalman filter state, all biases are set to zero, and the calibration matrix is set to the identity matrix. Note that the estimated calibration parameters have no unit since the magnetometer output values are normalized by the manufacturer. All filters use 3,000 particles and are updated at a rate of 10 Hz. The average processing time of all filters is well below the length of the data set. Therefore, the implementation is capable to run in real time.
3.4 Attitude Errors
During the derivation of the filter equations in Section 2.4, the assumption was made that the train frames of the Diesel railcar and the steam locomotive coincide 
and that the calibration parameters can be modeled independently from the attitude of the trains. To show that this is the case for the presented measurements, the attitudes for both trains were estimated with a loosely coupled INS/GNSS algorithm. In particular, the measurements of a KVH 1750 IMU were combined in an error-state Kalman filter with the position and velocity measurements of a u-blox LEA-M8T GNSS receiver. Figure 4 shows the difference in the roll and the pitch angles in degrees over the along-track position. The shown differences are the Euler angles derived from the rotation matrix, 
, between the two frames. In the figure, it is clearly visible that both trains have a comparable attitude along the whole track with errors below 0.5° for most parts.
Difference between the roll and pitch angles between the railcar and the steam locomotive
The potential error the misalignment introduces into the measurements is compared to the measured magnetic variations of the steam engine in Figure 5. The error is the value that one would measure when considering a sensor axis to be orthogonal to the measured magnetic vector while, in reality, it is misaligned. For a magnetic magnitude of 50 μT, the error due to a roll angle difference, therefore, is ±50 μT sin(ΔRoll). The magnetic variations are obtained by subtracting the moving average over the last 50 m of the corresponding axis from the magnetic map. The comparison in Figure 5 shows that, compared to the error, the variations are significantly larger for many parts of the track, leading us to the conclusion that neglecting the misalignment seems reasonable.
Magnetic variations (blue) recorded with the steam engine after removal of the mean and the attitude error ±50 μT sin(ΔRoll) (red)
3.5 Results
In this section, the results of the position and parameter estimation are presented. For better comparability, the full and reduced SLAC approaches are compared to a particle filter with fixed calibration parameters and a particle filter that uses no calibration at all. For the filter with fixed parameters, first, the magnetic field values from the map are obtained with the GNSS position and, second, the parameters are estimated with a least-squares estimator. The estimated parameters are then considered constant and known to the particle filter. The calibration is performed with the magnetometer measurements the train recorded on the first 500 m of the track. For the reduced SLAC algorithm and the particle filter without calibration, the magnetometer data of the railcar was only rotated with Equation (24) into the same body frame in which the map was also recorded. The full SLAC approach directly uses the raw magnetometer measurement without prior rotation. Due to the probabilistic nature of particle filters, all filters are evaluated over 100 Monte-Carlo runs.
3.5.1 Position Error
For error calculation, the minimum mean square error estimate (MMSEE), 
, of the along-track position:
28
is calculated from the posterior Dirac mixture density where 
is the along-track position of particle i. The ground truth of the along-track position sk is the value for which the distance between the ECEF position associated with it in the map and the GNSS ECEF position at time step k is minimal. Hence, the along-track error is simply the difference between the ground truth and that MMSEE, 
. The resulting along-track error of the four different filters is shown in Figure 6.
Along-track position error for the four different filters; for the pre-calibrated and the pre-rotated particle filter, the errors are shown only until the time instant at which the error reaches a magnitude of 30 m.
The shown errors for the two SLAC approaches are those resulting from the Monte-Carlo run with the highest root-mean-square error (RMSE), although the error showed similar behavior for all runs. For the remaining two filters, the error for the run with the smallest RMSE was chosen. From Figure 6, it becomes clear that, even for the best runs, both filters diverge and the position error increases up to kilometers. To not clutter the figure, the error is shown only until the time it reaches a magnitude of 30 m for the first time. Beyond that, the error is no longer shown. The particle filter without calibration immediately diverges after the start and is not able to track the train position. The pre-calibrated filter is able to perform localization for ≈ 260 – 275 s, depending on the run, but then loses track. At the beginning, the filter position error is smaller than for SLAC due to its known calibration parameters. After the parameters are estimated by SLAC and before the pre-calibrated filter diverges, the errors are almost identical. This shows that, after an initial phase with larger errors, the performance of SLAC is close to a filter that uses a calibrated sensor. While we expected that the particle filter with uncalibrated data would diverge, the divergence of the pre-calibrated filter was not expected. We suppose that this was caused by changes in the parameters and violations of the linear assumption for the calibration model. This will be discussed in more detail in the next section.
The SLAC algorithms have a comparable performance, but overall, the reduced SLAC shows a faster convergence in the beginning, also reducing the magnitude of the error. This is also reflected in the RMSE. For SLAC, the average RMSE over all runs is 2.44 m and, for reduced SLAC, is 1.85 m. In comparison, the two particle filters without simultaneous calibration diverged in every run and achieved RMSE values above one kilometer. At around 330-350 s in Figure 6, the error grows up to ≈ 30 m and is then reduced again to a few meters. The larger error occurs at a part of the track where only small variations in the magnetic field are observed as can be seen from the magnetic data in Figure 7. Since the filter uses only magnetometer measurements and a motion model for localization, the lack of sufficient or unique magnetic variations is problematic (i.e., in such scenarios, the filter could diverge). For applications that require high reliability, it is therefore strongly recommended to consider an odometer to aid the estimation. Due to the generic nature of the particle filter, incorporating additional sensors is straightforward.
Measured magnetic field variations: The railcar measurements are rotated to the body frame of the steam engine. The shown data is the output of the Xsens sensor which is normalized to 40 μT.
3.5.2 Calibration Parameters
The calibration parameters estimated by the full and reduced SLAC are depicted in Figure 8 and Figure 9. The estimated parameters applied to the map data are shown in Figure 10 and Figure 11. Note, in contrast to Figure 10, in Figure 11, the rotated railcar magnetometer data is plotted, which is also the input to the reduced SLAC algorithm. From Figure 10 and Figure 11, it can be seen that, after an initial phase of larger errors, both algorithms are able to closely fit the map data to the measured data.
Estimated parameters by the full SLAC approach
Estimated parameters by the reduced SLAC approach
First 200 s of measurement data and calibrated map data for the SLAC algorithms; the calibrated map data is obtained with the MMSEE of the train position and calibration parameters.
First 200 s of measurement data and calibrated map data for the reduced SLAC algorithms; the calibrated map data is obtained with the MMSEE of the train position and calibration parameters.
To our surprise, in Figure 8 and Figure 9, instead of converging to a constant value, the calibration parameters kept changing while the railcar drove along the track. This is in contrast to what we saw from earlier experiments with a model train presented in Siebler et al. (2021b) where the parameters converged towards a constant value. For the majority of the time, the parameters were slowly drifting but also a few more sudden changes were observed. Unfortunately, it was not possible to identify a single cause for this behavior. As for the model train experiments in Siebler et al. (2021b), the same sensor type was used, so the changes observed in Figure 8 and Figure 9 most likely were not caused by internal sensor parameters. We suspect that the changes were caused by a combination of violated model assumptions, unavoidable imperfections, and possible changes in the measurement setup. For the linear model assumption to hold, the surroundings of the sensor that influences the measurements must be completely static (e.g., no moving ferromagnetic material or power lines with slowly changing DC currents). This is hardly fulfilled for a train due to its steel wheels, gearboxes, and electrical systems. In addition, as the sensors in the steam locomotive and the railcar were not mounted at the exact same height and did not have the same exact roll and pitch angles, this could result in a nonlinear or time-variable relation between both sensor measurements. Another possible source of error is the map. The map most likely contains errors caused by the GNSS measurements providing the ground truth during map creation. For example, when the magnetic field in the map has a sudden shift compared to the true magnetic field, the calibration parameters are adapted by the Kalman filter to partially compensate for this error at the moment a particle reaches this position. Comparing Figure 8 and Figure 9 with the data in Figure 7, we see that the sudden parameter changes occur when the measurements have strong variations. Between two strong variations, the parameters are drifting, indicating that the parameter estimation is not accurately possible due to a low signal-to-noise ratio. When the particle cloud reaches a position with strong variations, smaller errors in the parameters lead to bigger residuals in the Kalman filter and the parameters are quickly corrected to reduce the residuals again.
The combination of model mismatch and the errors in the map could explain why a constant set of calibration is not sufficient to accurately calibrate the magnetometer for the whole data set and why the pre-calibrated particle filter diverged in all Monte-Carlo runs.
4 CONCLUSION
For magnetic field-based localization, magnetometer calibration plays an essential role. Therefore, we proposed in Siebler et al. (2021a, 2021b) a simultaneous localization and calibration (SLAC) algorithm based on a Rao-Blackwellized particle filter that exploits the conditional linearity of the calibration parameters. This reduces the number of required particles and lowers the overall complexity. To further lower the complexity, in this paper, a reduced version of the SLAC algorithm was introduced that only estimates a subset of parameters. In the reduced SLAC algorithm, the dimension of the Kalman filter was reduced and only one scale factor and bias per sensor axis was estimated. The old and new SLAC algorithms were evaluated with a data set recorded with a Diesel-powered railcar and a magnetic map recorded with a steam locomotive. Overall, an average RMSE position error below 2.5 m could be achieved for both SLAC versions. Due to faster convergence, the reduced SLAC algorithm achieved a lower RMSE, but the overall performance of both algorithms was similar. The results clearly show that SLAC also enables magnetic localization for trains that use uncalibrated sensors and are different to the train that was used to create the magnetic map.
HOW TO CITE THIS ARTICLE
Siebler, B., Lehner, A., Sand, S., & Hanebeck, U. D. (2023). Simultaneous localization and calibration (SLAC) methods for a train-mounted magnetometer. NAVIGATION, 70(1). https://doi.org/10.33012/navi.557
ACKNOWLEDGMENTS
We thank the HSB for their support during the measurements.
This work is partially carried out in the Ubiquitous Spatio-Temporal Learning for Future Mobility (ULearn4Mobility) project funded by the Helmholtz AI Cooperation Unit under the grant ZT-I-PF-5-49 and the Intelligent Magnetic Positioning for Avoiding Collisions of Trains (IMPACT) project funded by the Bavarian Ministry of Economic Affairs, Regional Development and Energy under the grant DIK-2002-0016 / DIK0175/01.
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.
Jump to section
Related Articles
Cited By...
- No citing articles found.