Graph-based fusion of magnetometer and learned-based inertial odometry

—The 3D position estimation of pedestrians is a vital problem in the development of virtual reality, augmented reality, and the internet of things. The learning-based inertial odometry is a very potential auxiliary method of pedestrian positioning due to its low position drift and immunity to external environmental inﬂuences. However, in many cases, the drift error of the heading is still the main factor that causes the rapid divergence of the position estimated by the learning-based inertial odometry. This paper proposed a graph-optimization-based estimation method to fusing learned-based inertial odometry and magnetometer measurements for obtaining lower drift position. The proposed algorithm does not need to calibrate the magnetometer bias, and effectively resist the inﬂuence of magnetic interference in the indoor environment, and can provide a very reliable absolute magnetic heading correction. Test results show that the proposed method can obtain better positioning performance than other methods using calibrated magnetometer observations.


I. INTRODUCTION
Pedestrian positioning plays a vital role in the internet of things (IoT), location-based services (LBS), and augmented reality (AR) [1]. Exiting high accuracy positioning techniques generally rely on LiDAR, camera, or wireless sensors. LiDARbased simultaneous localization and mapping (SLAM) [2] [3] which relies on professional equipment, can provide stable and reliable high-precision positioning. However, the LiDAR is heavy and power-hungry and is not suitable for pedestrians. Camera-based SLAM include visual odometry [4] and visualinertial odometry [5] [6] based built-in camera and inertial sensor of the smartphone is a low-cost and high-percision positioning solution with the potential to be widely used. Unfortunately, VIO still cannot work well in low-textured and dynamic illumination environments [7]. Wireless sensorbased methods (e.g., Ultra-wideband, 5G) rely on pre-arranged signal base stations. The cost is extremely expensive, which is unrealistic for a large-scale indoor environment.
Inertial Measurement Unit (IMU) can be used for providing 3D motion estimation independent of the external environments. Moreover, almost all wearable smart terminals (e.g., phones, watches, and glasses) are equipped with IMU. Therefore, from the perspective of system cost and sensor properties, the IMU-based method is ideal for providing pedestrian position in the actual application. This function is known as strapdown Inertial Navigation System (INS) [8]. However, the strapdown INS, which uses consumer-grade IMU, can not individually keep the accuracy of the position for more than several seconds. Because the position accumulation error of strapdown INS is proportional to the square of time and the errors (e.g., biases, scale factor, noise) of consumer-grade IMU are large and unstable.
Thus, pedestrian dead reckoning (PDR) is used as an alternative to providing usable position estimation of pedestrians. PDR uses prior knowledge of human motion patterns to mitigate accumulation error in PDR. More specifically, the PDR based on IMU usually uses the pedestrian gait motion pattern. Through identificate and classificate step, the estimated step length could be adopted to correct the estimated velocity of INS [9] [10] [11].
However, the motion pattern of pedestrians during walking is hard to be detected in real applications. Thus, the learnedbased PDR methods are proposed. This type of methods estimate velocity or displacement using raw IMU measurements, and show impressive accuracy compare to traditional PDR in experiments [12] [13] [14]. The learned-based PDR can be categorized into two groups.
The first groups of methods uses the magnetometer measurements, including IONet [12], RoNIN [13] and IDOL [15]. IONet estimates a 2D displacement using the IMU measurements represented in a global frame. The converting of IMU measurements from the device frame to the global frame uses the rotation estimated by the Android API. RoNIN is implemented based on a similar strategy but training and evaluating on a larger dataset. Meanwhile, more deep learning architectures include a temporal convolutional network (TCN), a residual network (ResNet), and long short-term memory (LSTM), which have been tested in field experiments. Both of the above two methods show better performance than traditional PDR. Therefore, the neural network inputs of IONet and RoNIN relied on absolute orientation estimated by Android API, and the Android estimate orientation uses magnetometer measurements. Thus, the position accuracy is easily degraded in areas with a highly perturbed magnetic field.
The second groups of methods are independent of magnetometer measurements, i.e., TLIO [14]. TLIO uses a multistate Kalman filter (MSCKF) to fuse the neural Network and the strapdown INS model. Unlike IONet and RoNIN, TLIO uses orientation estimated in the MSCKF to convert IMU measurements into a local gravity-aligned coordinate. This mechanism indicates that TLIO is independent of magnetometer measurements and not affect by magnetic perturbation. It also implied that the heading of TLIO is drift over time and affects the long-term positioning accuracy of TLIO. Although TLIO can provide accurate relative positioning results, the positioning error is still significantly affected by the drift of heading [15].
In theory, the magnetometer can provide absolute heading observation by measurement the earth's magnetic field. However, directly constraint heading uses magnetometer measurements is hard to work in indoor scenarios because the magnetometer measurements are not reliable [16]. More specifically, building materials and electronic equipment are highly perturbing the indoor magnetic field. Quasi-static magnetic field (QSF) technique only using magnetic constraint in the quasistatic field can effectively use magnetometer measurements to estimate heading in highly perturbed environments [17] [18]. The QSF can inhibit the increase of heading error, but the heading error is still slightly accumulated. This fact indicated that the heading drift still is the primary error source of PDR in long-term trajectory estimation. Furthermore, these methods can only use calibrated magnetometer measurements, which are hard to obtain in real applications.
The IDOL aims to solve this problem by training a neural network to estimate orientation based on inertial and magnetometer measurements. The IDOL shows the ability to capture magnetic field perturbation. Therefore, its' positioning accuracy outperformed the RoNIN and TLIO in the experiment. Meanwhile, the estimated orientation of IDOL can be utilized to improve the positioning accuracy of RoNIN. However, the learned-based orientation estimator may degrade in novel environments. Meanwhile, it is hard to collect data for every environment. On the other hand, this method can not work for magnetometer measurements with large biases.
To provide a global consistent orientation estimation, we hypothesize that the average magnetic field in a large area can reflect the absolute heading. In other words, we can provide a global consistent orientation estimation based on magnetometer measurements of a long trajectory. Furthermore, we model the magnetometer biases and estimate them online. Thus, the proposed method is robust to magnetometer biases. This paper has two major contributions: • We proposed a graph-optimization-based method to fuse learned-based inertial odometry (AI-IMU) and magnetometer measurements. This method can use magnetometer measurements in a long period to estimate absolute heading based on the average magnetic field in a large area. Thus, this method is robust to local magnetic field perturbation. • We model the magnetometer biases in the proposed method. So, the proposed method can work with an uncalibrated device, which is different from traditional methods. It is noticing that calibrate magnetometer manually before use is not user-friendly. We organize the remainder of the article as follows. Section II gives a brief description of the entire system. Section III gives a detailed description of the whole solution. Section IV uses field test to prove the proposed method shows improvement of accuracy. Section V summarized the entire study.

II. SYSTEM OVERVIEW
Compared with the traditional PDR, the learning-based inertial odometry shows gratifying positioning performance. However, the drift of the heading will still cause the performance of the learning-based inertial odometer to drop significantly. Then, the absolute heading angle based on the magnetometer observation is an available information source that can effectively control the heading angle error. As mentioned in the introduction, the use of magnetometer observations to train the network model cannot guarantee that the positioning performance will not decrease in the new building environment. The hypothesis that the main component of the environmental magnetic field is still the geomagnetic field is not limited to a single test environment. It can make the magnetic heading observation more universal. Therefore, this article combines the network model and the traditional empirical model to form a more robust pedestrian positioning solution. As shown in Fig. 1, this system consists of learned-based inertial odometry (denoted as AI-IMU) and a graph-optimization estimator.
The AI-IMU uses raw IMU measurement to provide 3D relative pose estimation, consisted of a neural network (denoted as Network in Fig. 1) and a Multi-State Cloning Kalman Filter (MSCKF). The neural network provides a 3D displacement between begin and end of a sliding window, and the sliding window length is 1 Second. The input of the neural network is IMU measurements in the sliding window represented in a local gravity-aligned coordinate frame. To obtain a continuous and stable pose estimation, MSCKF is adopted to fusing the information from IMU mechanization and the neural network (denoted as Network in Fig.1). INS mechanization is employed to predicts the system state (including pose, velocity, and sensor biases) and the corresponding uncertainty, which uses IMU raw measurements.
The graph-optimization estimator (denoted as Graph Optimization in Fig. 1) uses relative 3D pose coming from AI-IMU and magnetometer measurements to estimate trajectory without drift of heading. Considering the computational cost, the estimator only estimates keyframes and utilizes a sliding windows strategy. In detail, AI-IMU provides the relative pose and corresponding uncertainty between the current and previous keyframes, and the magnetometer provides magnetometer measurements of keyframes. The estimator builds up absolute and relative orientation constraints based on these magnetometer measurements to mitigate heading errors. Coordinate frames used in the proposed system. There are two groups of coordinate frames. The MSCKF uses F N , F B t and F L t . The graphoptimization estimator uses F G , F B t . Both F N and F G are gravity-aligned coordinate frames aligned with the center of IMU at the initial moment.
In this system, we implement a two-layer estimator by weighing accuracy and time consumption. As mentioned before, we aim to use the average magnetic field during a long period to provide an absolute heading constraint. To obtain an accurate estimation of absolute heading, we need to save states as more as possible in the method. However, the estimated system states in the system are limited due to providing motion estimation in real-time. Furthermore, the AI-IMU needs highfrequency corrected by neural network's output with highfrequency. Therefore, based on the high-frequency output of AI-IMu, we use a graph-optimization-based estimator to fuse AI-IMU results and magnetometer measurements at a low frequency.

A. Coordinate Definition
This paper defines two groups of coordinate frames shows in Fig.2: coordinates adopted in the MSCKF and coordinates adopted in the graph-optimization estimator. It is noticing that these two groups of coordinate frames are adopted in two submodule (AI-IMU and graph optimization) individually.
The MSCKF uses three coordinate frames: the navigation frame denoted as F N , the t-th body frame denoted as F Bt , and the t-th local gravity-aligned coordinate denoted as F Lt . The body frame aligned the coordinates of IMU at t moments. F N is a gravity-aligned coordinate. It is aligned with the IMU center at the initial moment. F Lt is a gravity-aligned coordinate frame. It has the same position and heading to F Bt . The MSCKF estimates the 3D motion and the sensor bias of IMU. The 3D motion parameterized as position (t nbt ) of the t-th body frame, rotation (R nbt ) from F Bt to F N and velocity The graph-optimization estimator uses two coordinates: the global frame denoted as F G , and the t-th body frame denoted as F B . The global frame is the reference frame of the estimator. It is a gravity-aligned coordinate and locates at the position of IMU at 0-th moment. The magnetic reference field is represented in F G and denoted as B G . Two adjacent keyframes' pose in estimator denoted as {R gb t−L , t gb t−L } and {R gbt , t gbt }. The L defined the distance between keyframes. In the implementation, it sets to 1 second.
The raw IMU measurements at t-th moment are denoted as a Bt t and ω Bt t . The magnetometer measurements represented in F Bt at t-th moment are denoted as B Bt t . Furthermore, the IMU measurements represented in F N denoted as a N t and ω N t .
B. AI-IMU 1) System State Definition: The full system state at t-th moment is defined as: where η is cloned system states, and s t is the current system state. s t and η are defined as R nbt is rotation from F Bt to F N . t nbt and v nbt are position and velocity of F Bt in F N . b a and b g are the IMU accelemeter and gyroscopt biases. In this paper, the IMU noise model defined asâ Bt t = a Bt t + b a + n a (4) a bt t andĝ bt t are the measured values of acceleration and angular rate, respectively. a bt t and g bt t are the true values of accelerating and angular velocity, respectively. n a and n are random noise variables following a zero-centered Gaussian distribution.
The MSCKF is an error-state-based indirect Kalman filter. Thus the MSCKF estimates the error state, which indicates the difference between estimated and real values. The error state defined as Hence, the dimension of the system is 15 + 6m, where m is the number of cloned system states and 15 is the dimension of s t .
In detail, the error state is defined aŝ The covariance of the error state are defined as Σ δs,δη1 represent covariance between s and η 1 .
2) State Propagation: The filter propagates the system state using the IMU raw measurements based on IMU mechanization. In pedestrian positioning, the walking speed of the pedestrian is low. Thus the change of gravity orientation caused by IMU movement can be ignored. In detail, the gravity g N is assumed to be equal during positioning. Furthermore, the earth's spin is ignored. We assumed that the IMU only measures the rotation relates to the navigation frame F N . The IMU mechanization utilized in this approach is defined as 3) Measurement Update: The measurement update of the AI-IMU uses the 3D relative displacement estimated by the neural network. The input of the neural network is IMU measurements represented in a local gravity-aligned coordinate frame. To simplify the implementation, we send the IMU measurements represented in the navigation frame F N to the network. These IMU measures can be calculated based on the following equations.
The input of network denoted as {â t−100:t ,ω t−100:t }. It means the latest 100 IMU measurements. Since the IMU samples 100Hz, the input is IMU measures in the latest 1 second. The neural network inference process defined as This equation cannot be directly adopted in the measurement function because the displacement represented in the navigation frame imposed an absolute heading constraint. In other words, measurement update based on this equation causes absolute heading observable. Nevertheless, the absolute heading is never observable in this system. To solve this problem, we convert this displacement to the local gravityalignment frame F Lt as mentioned before. The conversion of displacement and corresponding covariance defined aŝ 3. An illustration of the factor graph for fusing AI-IMU and magnetometer. AI-IMU provides the relative magnetic factor and gravity constraint factor. The magnetic factor and relative magnetic factor are established based on magnetometer measurements. The magnetic norm factor constraints the norm of the local magnetic vector.
Here,R yawt is the rotation matrix represent heading component ofR nbt . More specifically, we can decomposesR nbt to three individual rotation matrix.
R yawt ,R pticht , andR rollt denote yaw, pitch and roll respectively. Hence, the measurement function is defined as: . Furthermore, we employed a χ 2 -test to avoid abnormal results provided by the network block. 4) Relative Pose and Uncertainty: In the proposed system, the AI-IMU function as an odometer which output relative poses for further fusion. In detail, we have to provide relative pose between {R nb t−L ,t nb t−L } and {R nbt ,t nbt }. Covariance matrix of the relative pose is also necessary. The relative pose is defined as To modeling the uncertianty of this relative pose, we can calculate the corresponding covariance based on follow equation.
C. Graph-optimization Heading Fusion 1) Problem Definition: The graph-optimization estimator fuses relative pose and magnetometer measurements to provide trajectory without heading drift. The definition of the factor graph shows in Fig.3. Three states are estimated in the factor graph, including keyframe pose {R gbt ,t gbt }, magnetometer bias b B , and magnetic field B G . The magnetometer bias is the bias of the magnetometer. Through estimate magnetometer bias online, the estimator can work for uncalibrated magnetometers. Furthermore, we estimate the average magnetic field at the initialization stage and recognize this average magnetic field not changing over the whole experiment area. The initialization stage is described in Section III-C3. Otherwise, we only keep system states in a short period in the estimator. In detail, we maintain 300 latest states in the estimator and select one keyframe per second. So, the estimator contained state in the latest 300 seconds to estimate an average absolute heading. As mentioned before, we hypothesis that the average magnetic field of a large area can reflect the absolute heading. More specifically, the absolute heading estimated by the graphoptimization is according to the average magnetic field of the latest 300 seconds. Thus, this absolute heading should be close to the real absolute heading and not effect by local magnetic field perturbation. Furthermore, to eliminate the oldest keyframe without loss of information, we adopt the marginalization technique, which has been widely employed in the sliding window graph-optimization problem. Section III-C4 gives details of the marginalization technique.
The full state vector at t moment in the sliding window is defined as We estimate the maximum posterior estimation of the sliding window as the result of the estimator. In detail, we estimate X t by minimize the following cost function Σ represents the Mahalanobia norm of residual. ρ(·) is the Huber norm [19], defined as r prior represents prior constraint. r B G represents the whole strength constraint of B G . r g and r odo represent gravity orientation constraint and relative pose constraint, are provided by AI-IMU. r B and r ∆B represent absolute and relative magnetic measurement constraint, are based on magnetometer measurements. Details of residuals are provided in Section III-C2. The cost function (28) defines a non-linear least squares problem. We uses Levenberg-Marquardt algorithm [20] [21] which implemented in Ceres Solver [22] to solve this problem. More specifically, we linearize the (28) and solve the linearized equation iteratively.
2) Residual Definition: r prior represents prior constraint and defined as A and b are defined based on prior distribution. In detail, x can be recognized to follows normal distribution N (b, (A T A) −1 ). r B G represent the total strength constraint to the magnetic field. It is defined as Z B G represents norm of the local average magnetic field vector. In partical, this value can be manually selected based on theory model or average magnetic field estimated before. r odo represents residual of the relative pose. It is defined as t btb t−L and R btb t−L are relative position and rotation represented in F B t−L and calculated based on (25) and (24). Log SO(3) (·) represents logarithm function for SO(3) and outputs a three-dimensional vector. Furthermore, the covari- It is worth to be noticing that all linearization operations for rotation use the same formula. In detail, we use perturbation at the right side of the rotation matrix. r g represent constraint of gravity orientation. It is defined as g Bt and g G represent gravity vector in F Bt and F G respectively. g G is predefined since the global frame are gravity aligned coordinate frame. g Bt is calculated based on rotation estimated in AI-IMU. More specifically, it calculated based on Since the relative pose residual r odo uses relative rotation and displacement to constraint relative poses, the pitch and rolling angle relate to the navigation frame implicated in the MSCKF of AI-IMU are not maintained. Thus, we use the orientation of the gravity vector to constraint absolute rolling and pitch angle. r B represents absolute heading constraint based on magnetometer measurements. It is defined as B Bt represent magnetometer measurements at t moment. This residual provide constraint for absolute heading. Otherwise, the Huber norm is applied for r B to avoid effect of abnormal magnetometer measurements. The abnormal magnetometer measurements can caused by abnormal local magnetic field or other error sources. r ∆B represents relative heading constraint based on magnetometer measurements. It is defined as This residual provides the relative rotation constraint based on magnetometer measurements for adjacent keyframes. To avoid perturbation of the magnetic field, we also applied the Huber norm for this residual. Through combination of r ∆B , r B , and r B G , we can estimate b B and B G simultaneously. More specifically, b B and B G are observable when the IMU is rotated around any axis.

3) Initialization:
We aim to design a method that can maintain absolute heading related to the first short period. As we hypothesis that the average magnetic field in the long period can reflect absolute heading, the average magnetic field in the first sliding windows is equal to that in any other sliding window. So, we can estimate the average magnetic field for the first sliding window and use this average magnetic field to provide the absolute heading constraint in the rest of the trajectory.
In the initialization stage, all system states in X t are estimated. The absolute position and absolute heading still are unobservable in the system. To solve this problem, we fix the pose of the first keyframe ({R gb0 , t gb0 }) and set a prior constraint to this pose. Furthermore, we first optimize the (28) when there are more than ten keyframes contained in the sliding window. This strategy is helpful to avoid estimator estimate an abnormal magnetometer bias b B and magnetic field vector B G at the first time.
When the sliding window has been grown to a particular length, we fix the magnetic field vector B G . In this implementation, we fix B G when the sliding window contains 300 keyframes. In theory, through adopting the marginalization technique, we can get the optimal posterior estimation of the magnetic field vector B G . Nevertheless, the marginalization operation uses a Gaussian distribution to approximate the information of marginalized residuals. This strategy causes a loss of information. In this problem, the loss of information causes the magnetic field vector B G can be slowly changing over the whole positioning process. Thus, we can not maintain an absolute consistent heading if we do not fix the magnetic field vector B G . 4) Marginalization: To eliminate computation cost, we employ a sliding window estimator. We adopt the marginalization technique [23] [5] to give an approximate of information contained in removed residuals. It is widely used in visual odometry. When we aim to marginalize out a set of states δX m , we have to adopt the Schur-Complement operation to the sub-problem contained by X m and relative system state X r . In detail, we linearize the sub-problem at particular linearization point to the linear system defined as X op m and X op r represent linearization point of X m and X r . In detail, (37) is the linear system used in Gaussian-Newton method. Through applicate the Schur-Complement operation to (37), we get the linear system not contained X m .
(40) can be easily convert to a prior constraint for remained system states X m . In the view of probability, X m − X op m follows a zero-mean Gaussian distribution. Fig. 4 illustrates the marginalization strategy in the proposed method. In our approach, we eliminate the pose of the oldest keyframe. Relative system states contained magnetic field vector B G , magnetometer bias b B , and pose of keyframe next to the oldest one. The relative residuals contained r prior , r odo , r g , r B and r ∆B . After eliminating the oldest keyframe, the prior distribution constraint for the relative system states is added to the problem.

IV. EXPERIMENTS
In this section, we show the magnetic field of the experiment scenarios. Then, we compared the proposed method to other magnetic fusion methods and verified that the proposed method could work with uncalibrated magnetometer measurements.
The remainder of this section is organized as follows. Section IV-A describs the test implmentation details. Section IV-B describes the metrics to evaluate the performance of the proposed method. Section IV-C shows magnetic fields and compares the accuracies of all methods.

A. Experiment Description
An Asus Tango phone is mounted to the head for data collection. The Tango phone contains a global shutter fisheye camera, a depth camera, an IMU, and a magnetometer. And the IMU measurements, magnetometer measurements, and ground-truth poses are collected at 100 Hz. The magnetometer measurements are the calibrated measurements provided by Android API. Since the change of the phone's working state will form an additional equivalent magnetometer bias, we use the ellipsoid fitting method to manually compensate the magnetometer bias to obtain a more accurate magnetic field vector. The ground truth poses are outputs of visual-inertial odometry. Indeed, the tango phone can use the embedded camera and IMU to estimate 3D motion through the visualinertial odometry technique. Furthermore, it can build a visual map for the selected area and do positioning in this map. So, we use the result of positioning in a pre-built map as ground truth if possible. Otherwise, in the senario the tango phone can not build a map, we select trajectories that visualinertial odometry can provide good positioning accuracy. We use outputs of visual-inertial odometry at ground truth directly.
We use the AI-IMU method described in LLIO [24] which is a lightweight version of TLIO [14]. We use the ResMLP256 LLIO-Net [24] to estimate short period displacement in the AI-IMU. The LLIO-Net is trained used 40 hours dataset. Six persons capture the dataset use three devices in two buildings. It is worth noticing that the magnetometer measurements were never used during training the LLIO-Net.
We compare the proposed method with several other methods. Details of these methods are given as follows.
• AI-IMU The AI-IMU, as mentioned in Section III-B, is pure learned-based inertial odometry without dependency on magnetometer measurements. It is the basis for all the methods described below.
The Mag Filter is the modification of AI-IMU, adding a measurement function used magnetometer measurements. In detail, it measurement update defined as B N is magnetic field vector in the navigation frame F N . It is directy given in the experiments. B Bt is magnetometer measurements at t moment. An χ 2 -test are employed to avoid abnormal magnetometer measurements.
The QSF Filter uses magnetometer measurements based on quasi-static field (QSF) hypothesis. It detect quasistatic magnetic field, and adopting magnetic constraint in the quasi-static magnetic field. More specifically, it detect quasi-static magnetic field based on detecting changing of magnetic strength [17]. In a quasi-static magnetic field, its measurement update defined as Here, B QSF is magnetic field vector of current quasistatic magnetic field. It is defined as k moment is the first moment of the current quasi-static magnetic field.
The Graph-Based is the proposed graph-optimizationbased fusion algorithm. In this approach, we select one keyframe per second and keep a sliding window contained 300 keyframes. We denote the ground-truth trajectory as GT in the rest of the paper.
Since the Mag Filter and the QSF Filter do not model the bias of the magnetometer, we give the calibrated magnetometer measurements. Thus, these two methods are only affected by the perturbation of the environmental magnetic field. Furthermore, we manually added a bias to the calibrated magnetometer measurements and adopted it in the Graph-Based method. This strategy is adopted to verify that the Graph-Based method is robust to magnetometer bias. At this condition, the Graph-Based method needs to estimate absolute heading and magnetometer bias simultaneously.

B. Evaluation Metric
In this paper, we use aligned average position error in the horizontal plane to evaluate the positioning performance of the above algorithms. Position of estimated trajecotry and ground truth trajectory are represented in navigation frame are denoted as {t 0 , · · · ,t N } and {t 0 , · · · , t N }.R t and R are 2D rotation matrix.t t and t are 2D translation. We use the first 20% of the estimated trajectory to align to the ground truth trajectory. In detail, we estimate the 2D rotation and translation {R a , t a } fulfill the below formular.
Then, the average position error is defined as R a and t a are alignment rotation and translation described in (46). · represents l2-norm.

C. System Performance
We described a systematic comparison between the proposed system and other methods in indoor positioning scenarios. As mentioned in Table I, we use six trajectories at different places and with different shapes to evaluate the magnetic fusion result. The total length of trajectories except trajectory A is around 700 meters. All trajectories except trajectory A walk more than 10 minutes to make the drift of heading causes significant positioning error. Trajectory A is walking in a small office building where the Tango phone can construct a prebuilt map. The position error is contained in Table I. The proposed method (Graph-Based) shows the minor positioning error over all cases. The position error of Graph-Based (the proposed method) is significantly smaller than other methods. Fig. IV-C shows the cumulative distribution funtion (CDF) of position error of all trajectories. The performance of Graphbased is better than others. Furthermore, the performance of Mag Filter and QSF Filter is not stable. The QSF Filter shows better accuracy than the Mag Filter in Trajectory A, B, C, and D. In trajectory E and D, Mag Filter shows better accuracy than QSF Filter. Fig. 6 and Fig. 7 give detail of the magnetic field and comparision between algorithms. Fig. 6 (a) and Fig. 7 (a) show the local magnetic field vector. It is calculated based  on the magnetometer measurements B Bt and the rotation matrix R nbt of ground truth trajectory. The local magnetic field B N local defined as The magnetic field is roughly consistent but perturbated at some positions. Thus, the hypothesis that the average magnetic field in sub-areas can reflect the absolute heading is fulfilled. Fig. 6 (b) and Fig. 7 (b) show angle difference between each method and ground truth. We calculate the angle difference based on the raw outputs of each method. In detail, we provide initial rotation and position based on ground truth. Thus, the initial heading error is zero. To illustrate the effect of the local magnetic field heading, we provide the heading difference between the local magnetic field vector to the reference magnetic field vector. This heading difference is denoted as magnetic heading error in Fig. 6 (b) and Fig. 7 (b).
Since Mag Filter directly uses magnetometer measurements in the measurement update, its heading is easily be affected by the magnetic field perturbation. Meanwhile, its heading is slightly better than the heading of the magnetic field heading. The Kalman filter estimates the current heading based on the compressed prior information and current magnetic field heading simultaneously. In other words, its heading is the weighted average of previous magnetic headings. As shown in Fig. 6 (c) and Fig. 7 (c), its position accuracy is significantly degraded because of the heading perturbation.
QSF Filter is robust to magnetic field perturbation since it does not use magnetic constraint when the magnetic field variant. The heading error of the QSF Filter is increasing over time because the QSF does not use a global consistent reference magnetic field vector. Although QSF is not easily affected by local magnetic field perturbation, its heading error is higher than Mag Filter in the long term.
The proposed method (Graph-Based) is robust to magnetic field perturbation and does not drift over time. It shows the best relative accuracy of heading estimation. However, it may keep a bias of heading related to the heading of ground truth. This bias is because of the difference between the global average magnetic field and the average magnetic field at the beginning moment. In other words, the average magnetic field during the first sliding window contains a bias to the other areas. We evaluate the position accuracy use the aligned trajectory to avoid this problem. As shown in Fig. 6 (c) and Fig. 7 (c), the aligned trajectories show that Graph-Based provides the best relative positioning accuracy in the experiments. In the experiments, Graph-based, which uses uncalibrated magnetometer measurements, shows better performance than other methods which use calibrated magnetometer measurements.
Furthermore, the proposed method can do positioning even the magnetometer measurements with a bias. As mentioned before, Graph-Based uses magnetometer measurements with a known bias in the test. The bias is [30, −20, 30]µT . Table  I gives the estimated magnetometer bias of Graph-Based at each trajectory. Fig. 6 (d) and Fig. 7 (d) give the variation of estimated magnetometer bias over time in trajectory A and B. We estimate the y and z components of the magnetometer bias close to the given magnetometer bias. The x component is convergent slowly or even not convergent to the given value. That phenomenon may be caused by the x-axis roughly oriented to up. During the walking process, the pitch and rolling angle are merely variants. Thus the observability of the x-axis is weak. In some cases, the x component of magnetometer bias may be unobservable without the constraint of magnetic field vector norm B G . The x component does not significantly influence the heading estimation when its observability is weak. Table II gives the estimated magnetometer bias and position error of Graph-Based at trajectory A, B, and C. The first column shows the magnetic biases in magnetometer measurements given to Graph-Based optimization. Graph-Based shows similar performance at different given magnetic biases. This fact indicates that the proposed method is robust to magnetometer bias.
Since the proposed method aims to run on mobile devices, its computational cost should be considering. Table III shows the time consumption of each block of the proposed method (Graph-Based) during run trajectory F. The total time length of trajectory F is 866 seconds. The proposed method costs 32.38 seconds on a computer with Intel i7-10700K. It runs at 27× of real-time on the computer. The proposed method combines AI-IMU with an additional graph-optimization-based estimator to fusing magnetometer measurements. AI-IMU consists by propagation (Section III-B2), measurement udpate(Section III-B3) and network inference ((18)) blocks. The bottleneck of AI-IMU is the network inference module which costs 18 seconds. The grah-optimization-based estimator denoted as Graph-Optimization in Table III. This block, which is described in Section III-C, contains a build graph, optimize the graph, and marginalization blocks. It costs 11.5 ms per time but only runs at 1Hz. Thus its total cost time is 30 seconds which is least than AI-IMU blocks. In summary, the added graph-based optimization block's time consumption does not significantly exceed that of AI-IMU. Meanwhile, the proposed method's computation resource-demanding is acceptable for running on mobile devices.

D. Discussion
The robustness for magnetic field perturbation is the main advantage of the proposed method. As shown in Fig. 6 (a) and Fig. 7 (a), the hypothesis that the average magnetic field in a large area can reflect the absolute heading is satisfied in these two experiment environments. Indeed, this hypothesis is satisfied in all 6 experiment trajectories collected in 4 different buildings. Under this condition, the proposed method can constrain absolute heading in areas with local magnetic field perturbation and provide better position accuracy than others, as mentioned in Section IV-C. However, since the sliding window length is limited, this hypothesis may not be satisfied in some cases. For example, when the pedestrian walks around a small area with similar magnetic field biases, these biases cannot be mitigated by estimating the average magnetic field of the sliding window. Therefore, those exceptional cases may not frequently be present.
Another advantage of this method is that the magnetometer bias is estimated online, as mentioned in Table II. Thus, this method can use for magnetometers, which are not calibrated or not well-calibrated, such as the one installed in smartphones or AR headsets. Indeed, this situation is common on consumergrade devices.
The proposed method is a combination of the AI-IMU module and the graph-optimization module. On the one hand, this design helps the proposed method achieve a trade-off between accuracy and computational cost. The AI-IMU module output high-frequency poses estimation, and the graphoptimization module run at a low frequency to fuse outputted relative poses of AI-IMU and magnetometer measures. Thus, the proposed method achieved significantly positioning accuracy improvement than AI-IMU but only slightly increased computation cost, as mentioned in Table III. On the other hand, this combination balances dependence on the dataset size and positioning performance. The AI-IMU can provide better relative position estimation accuracy than traditional PDR methods. It only uses IMU measurements which indicates that its performance only depends on the coverage of motion patterns in the training dataset and not on the coverage of the magnetic field perturbation in the dataset. The Graphoptimization module fuses magnetometer measurements by physical model rather than data-driven method. In summary, the proposed method could work in novel environments contain novel magnetic field perturbation. Meanwhile, it is easier to make dataset contains enough motion pattern than contains all magnetic field perturbation.
V. CONCLUSION In this paper, we propose a graph-optimization-based estimator to fuse AI-IMU and a magnetometer. The proposed method can use the absolute heading based on magnetometer observations to improve the pose estimation performance of AI-IMU in indoor environments. The reasons include: 1) the proposed method assumes that the mean value of magnetic interference is zero and uses a sliding window to suppress the effect of magnetic field perturbation. 2) the proposed method estimate magnetometer biases online to deal with the magnetometer bias change caused by the working state of the smartphone. The test results show that the proposed method, even if the uncalibrated magnetometer observations are used, still achieves the best position accuracy in all scenarios than the existing methods. Since the typical problems of the magnetometer (e.g., indoor environmental magnetic interference and frequent changes of the magnetometer bias) have been resolved, the proposed method can transform the magnetometer from an unusable state to a very useful state. Moreover, to meet the requirement of real-time positioning, we analyze the time consumption of each module of the proposed method. The test results show the proposed method can run 27× faster than in real-time and will not significantly exceed the time consumption of AI-IMU.
In future, we will focus on the tightly coupled fusion of the neural network, IMU mechanism, and magnetometer measurements for obtaining a more stable and higher-performance pose estimation. Moreover, learned-based magnetic field perturbation detection is also worth to be explored. ACKNOWLEDGMENT