Use of Magnetic Field for Mitigating Gyroscope Errors for Indoor Pedestrian Positioning †

In the field of indoor pedestrian positioning, the improved Quasi-Static magnetic Field (iQSF) method has been proposed to estimate gyroscope biases in magnetically perturbed environments. However, this method is only effective when a person walks along straight-line paths. For other curved or more complex path patterns, the iQSF method would fail to detect the quasi-static magnetic field. To address this issue, a novel approach is developed for quasi-static magnetic field detection in foot-mounted Inertial Navigation System. The proposed method detects the quasi-static magnetic field using the rate of change in differences between the magnetically derived heading and the heading derived from gyroscope. In addition, to eliminate the distortions caused by system platforms and shoes, a magnetometer calibration method is developed and the calibration is transformed from three-dimensional to two-dimensional coordinate according to the motion model of a pedestrian. The experimental results demonstrate that the proposed method can provide superior performance in suppressing the heading errors with the comparison to iQSF method.


Introduction
High accuracy and infrastructure-free indoor pedestrian positioning is a highly desired ability for finding and rescuing firefighters or other emergency first responders, or for personal navigation assistance. With the adoption of Global Navigation Satellite System (GNSS), positioning in outdoor environment has witnessed a dramatic success. However, GNSS performs poorly in indoor environment and is not suitable for indoor positioning because of weak signal reception and radio wave multi-path effect. Developers have suggested various indoor positioning technologies such as Ultrasound [1,2], Radio frequency (RF) [3,4] and Ultra wide band (UWB) [5,6]. The Ultrasonic positioning system localize a pedestrian based on the propagation time of ultrasonic wave and it can achieve an accuracy within a meter. However, it needs large-scale basic infrastructures resulting in extra costs. The RF positioning technology is one of the most cost-effective methods because wireless local area network (WLAN) is available in most indoor scenarios. Nevertheless, due to the multipath arising from reflection and the scattering characteristic of indoor environments, RF do not achieve high levels of accuracy when applied to large areas [7]. The accuracy of most of the RF based positioning systems is typically in the range of a few meters. Wi-Fi fingerprinting positioning is another solution for indoor positioning. In this method, first of all, a radio map will be created by collecting fingerprints of a scene and then it will be used to estimate user location [8]. The radio map has to be acquired by system platform and shoes. The calibration of the magnetometer is carried out before it is adopted for positioning.

Heading and Inclination Angle Estimation Using Magnetic Field
The Earth's magnetic field can be precisely measured and utilized to estimate heading by magnetometers. The three-dimensional Earth's magnetic field vector is shown in Figure 1.
where g denotes acceleration of gravity. Noting that the heading estimation proceeds only during the stance phase, when k r and k p in Equation (3) can describe the attitude of the sensors. During Here B denotes the total magnetic field, H is the horizontal component of B, D is the declination angle between H and the true north. B x , B y and B z are the three orthogonal components of total magnetic field. The magnetically derived heading with respect to the true North is estimated as: Inclination angle I is given by From Equation (1), it can be observed that the magnetic field components should be transformed to navigation frame (defined to be North-East-Down on the ground) so as to find the horizontal component of the Earth's magnetic field. Let (m k x , m k y , m k z ) be the sample of the three-axis magnetometer's measurement at time k in sensor body coordinate frame. Let (r k , p k , y k ) be the Euler rotation angles between body coordinate and navigation frame, then the transformation can be expressed by roll (r k ) and pitch (p k ) angles: where r k and p k are obtained from the accelerometer measurement as follows: where g denotes acceleration of gravity. Noting that the heading estimation proceeds only during the stance phase, when r k and p k in Equation (3) can describe the attitude of the sensors. During the stance phase, the heading y k is the only Euler rotation angle unknown. In addition, the magnetically derived heading ϕ is independent of the orthogonal field component B z according to Equation (1).
Therefore, three-axis magnetometer's measurement can be replaced by two-axis measurement and the replacement is given by In other words, the three-axis magnetometer can be used as the two-axis magnetometer in horizontal plane coordinate (North-East frame).

Magnetometer Error Source
Similar to inertial sensors such as accelerometers and gyroscopes, magnetometers also suffer from errors due to magnetic deviation which can result in an inaccurate heading measurement [25]. The mainly magnetic deviation can be divided into two categories: hard iron and soft iron distortions [26]. The presence of electro-magnetic systems and ferromagnetic materials beside the magnetometer are the main cause of these errors, that is the host platform itself may be responsible for them [27]. Hard iron distortions are caused by a magnetic source, which produce permanent bias added to each axis of magnetometer output. In other words, magnetic fields generated by different electronic sub-systems in the vicinity of the magnetometer are called hard iron magnetic sources. With the effects of hard iron distortions, the locus of magnetic field intensity will be a center-shift circle. Soft iron distortions are caused by much complex magnetic fields which are generated by ferromagnetic materials. The magnitudes of these magnetic field depend on the incident angle of the Earth's magnetic field on the material. Hence it changes as the system platform changes its attitude [28]. Soft iron distortions are non-linear and distort the magnetic field both in intensity and direction, making the locus of magnetic field intensity change from circle to ellipse.

Magnetometer Calibration Procedure
Many calibration algorithms have been proposed. One use measurements collected by a 360 • rotation of the levelled magnetometer in the horizontal plane to find the maximum and minimum values to estimate biases of the sensor. The iterative least squares algorithm is utilized to estimate the bias and the combined scale factors for a given data set [29]. This algorithm requires a proper uniform distribution of the data over all attitudes. The Constrained Least Squares Ellipsoid Fitting (CLSEF) method utilizes a constrained condition involved with the limited data to achieve an accurate ellipsoid fit [30]. Although these methods have been proved to be effective, they need large amount of computation. Furthermore, the performance of these methods will be degraded when applied to foot-mounted IMU platform, because the floor and shoes can cause interference to the Earth's magnetic field. To eliminate this interference, we propose an on-board magnetometer calibration approach.
When the foot is stationary on the floor, the platform is only affected by gravity. Then we can calculate the attitude angle of the platform according to Equation (3). Using the attitude angle, we can transform the magnetic field from sensor body coordinate system to navigation coordinate system. In addition, as abovementioned discussion, the orthogonal field component B z can be neglected during the heading estimation procedure, then the three-axis magnetometer is equivalent to the two-axis magnetometer in horizontal plane coordinate. The calibration procedure is described as follows:

Hard Iron Distortion Calibration Procedure
(1) Let the IMU be installed in the heel of the pedestrian's shoe, then the pedestrian walks around a loop as shown in Figure 2. Note that the rotation angle ψ should be as small as possible to store relatively complete measurements of the Earth's magnetic field. Here, we choose a reference value of ψ from 15 • to 35 • . (2) By using the stored data, the stance phase can be identified according to the Zero-Velocity detection algorithm in Reference [31].
(3) Transform the magnetometer's measurement from the sensor body coordinate frame to the navigation frame according to Equation (4). Then the average maximum and minimum values for each of the axes can be calculated. (4) Let max(x|y) and min(x|y) denote the maximum and minimum values of x and y axe respectively. Then, determine offsets of each axes as follows: offset_x = (max(x) + min(x))/2 offset_y = (max(y) + min(y))/2 (6) where offset_x and offset_y represent the offset of x and y axe respectively. Then these offsets are subtracted from the magnetometer's data that have been transformed at step (3) to eliminate the hard iron distortion effects.
store relatively complete measurements of the Earth's magnetic field. Here, we choose a reference value of ψ from 15° to 35°.
(2) By using the stored data, the stance phase can be identified according to the Zero-Velocity detection algorithm in Reference [31]. (3) Transform the magnetometer's measurement from the sensor body coordinate frame to the navigation frame according to Equation (4). Then the average maximum and minimum values for each of the axes can be calculated. (4) Let max( | ) x y and min( | ) x y denote the maximum and minimum values of x and y axe respectively. Then, determine offsets of each axes as follows: where offset _ x and offset _ y represent the offset of x and y axe respectively. Then these offsets are subtracted from the magnetometer's data that have been transformed at step (3) to eliminate the hard iron distortion effects.

Soft Iron Distortion Calibration Procedure
Compensation for the soft-iron distortion includes finding mathematical representation of the ellipse that is caused by disturbed magnetic field. The calibration procedure is as follows: (1) Find major axe and minor axe of the ellipse by a loop calculation of stored magnetic field database M : where ∈M ( , ) x y m m , max( ) r is the length of the major axe and min( ) r is length of the minor axe. For both max( ) r and min( ) r , we can get corresponding coordinates from M : x y m m .
(2) Determine inclination angle of the ellipse by: Then, the rotational matrix can be constructed to align the rotated ellipse with one of the coordinate system axes. The rotational matrix R is given by: cos sin sin cos (9) (3) Rotate the ellipse:

Soft Iron Distortion Calibration Procedure
Compensation for the soft-iron distortion includes finding mathematical representation of the ellipse that is caused by disturbed magnetic field. The calibration procedure is as follows: (1) Find major axe and minor axe of the ellipse by a loop calculation of stored magnetic field database M: where (m x , m y ) ∈ M, max(r) is the length of the major axe and min(r) is length of the minor axe. For both max(r) and min(r), we can get corresponding coordinates from M: (m x1 , m y1 ) and (m x2 , m y2 ). (2) Determine inclination angle of the ellipse by: Then, the rotational matrix can be constructed to align the rotated ellipse with one of the coordinate system axes. The rotational matrix R is given by: (3) Rotate the ellipse: where M rotated represents the set of coordinates after rotated. To scale the two-axis length of ellipse equal, the scale factor is required τ = min(r)/max(r). (4) Scale X axe (or Y axe) coordinate of the ellipse to make it circular: The coordinates set of the circle is rotated back to initial position using the rotational matrix and inclination angle: where M corrected represents the final result of the magnetometer calibration.

Quasi-Static Magnetic Field Detection
Although the Earth's magnetic field is a good measurement for heading estimation outdoor, it suffers severe perturbations in buildings. These perturbations induce random variations in the magnetic field estimations. As Figure 3 shows, the magnetic field magnitude in indoor environments is not constant due to the changing perturbations. Although the magnetic field indoor is not spatially constant, depending on the surroundings, it is possible to have locations as well as short periods when the perturbed magnetic field is constant which can be considered quasi-static [23]. The QSF method detects quasi-static magnetic field based on the change rate of total magnetic field magnitude. This method causes false-alarm detection easily when the magnitude of each axis changes with different ratio but the total magnetic field magnitude has very slight changes. The iQSF method adds two more detection constraints including the change rate of headings derived from magnetic field and angular rates (gyroscope output). The iQSF method can only take effect when the pedestrian walks along straight-line paths. If the pedestrian moves along non-ideal paths (curved or more complex paths), the change rate of the heading will be far greater than zero and the quasi-static magnetic field will be missing.
(4) Scale X axe (or Y axe) coordinate of the ellipse to make it circular: The coordinates set of the circle is rotated back to initial position using the rotational matrix and inclination angle: (12) where corrected M represents the final result of the magnetometer calibration.

Quasi-Static Magnetic Field Detection
Although the Earth's magnetic field is a good measurement for heading estimation outdoor, it suffers severe perturbations in buildings. These perturbations induce random variations in the magnetic field estimations. As Figure 3 shows, the magnetic field magnitude in indoor environments is not constant due to the changing perturbations. Although the magnetic field indoor is not spatially constant, depending on the surroundings, it is possible to have locations as well as short periods when the perturbed magnetic field is constant which can be considered quasi-static [23]. The QSF method detects quasi-static magnetic field based on the change rate of total magnetic field magnitude. This method causes false-alarm detection easily when the magnitude of each axis changes with different ratio but the total magnetic field magnitude has very slight changes. The iQSF method adds two more detection constraints including the change rate of headings derived from magnetic field and angular rates (gyroscope output). The iQSF method can only take effect when the pedestrian walks along straight-line paths. If the pedestrian moves along non-ideal paths (curved or more complex paths), the change rate of the heading will be far greater than zero and the quasi-static magnetic field will be missing.  To solve this problem, we use the rate of change in difference between the two headings mentioned above. As the error accumulation of the heading derived from gyroscope is very small during short periods, the heading derived from angular rate can be utilized to assess the stability of the magnetically derived heading. If the rate of change in differences between the magnetically derived heading and the heading derived from angular rate is approximate to zero, then the magnetic field can be considered as quasi-static. Herein, the detection of the quasi-static magnetic field is unrelated to the shape of the paths. In addition, the measurements of total magnetic field magnitude and inclination angle are also adopted. The magnetic field can be detected as quasi-static if both the To solve this problem, we use the rate of change in difference between the two headings mentioned above. As the error accumulation of the heading derived from gyroscope is very small during short periods, the heading derived from angular rate can be utilized to assess the stability of the magnetically derived heading. If the rate of change in differences between the magnetically derived heading and the heading derived from angular rate is approximate to zero, then the magnetic field can be considered as quasi-static. Herein, the detection of the quasi-static magnetic field is unrelated to the shape of the paths. In addition, the measurements of total magnetic field magnitude and inclination angle are also adopted. The magnetic field can be detected as quasi-static if both the changing rates of total magnetic field magnitude and inclination angle are equal to zero. The proposed detection approach is elaborated as follows.
Let θ m k be the magnetically derived heading at the k th step, θ g k represents that derived from the angular rate. B k denotes the magnitude of the magnetic field, I k denotes the inclination angle. The model of the measurement is as follows: where T and T denotes the transpose operation. Moreover, s θ k = θ m k − θ g k denotes the difference between the magnetically derived heading and the heading derived from angular rate at the k th step, s I k = |I k − I k−1 | denote change in the inclination angle of two adjacent steps; T denotes the measurement noise and we assume that the noise is independent identically distributed zero-mean, Gaussian with covariance matrix: where σ 2 θ , σ 2 I and σ 2 B are the noise variance of s θ k , s I k and s B k , respectively. {·} denotes the expectation operation. Mathematically, the magnetic field can be detected as a quasi-static field if s θ k = 0, s I k = 0 and s B k = 0. Then the detection problem can be formalized as a binary hypothesis testing problem and the detector can choose between the two hypotheses H 0 or H 1 defined as follows: where the hypothesis H 0 , H 1 denote non-static field and quasi-static field, respectively. Γ n denotes the step interval and Γ n = { ∈ N : n ≤ ≤ n + N − 1}. The performance of the detector is specified by the false-alarm probability P FA and detection probability P D . Neyman-Pearson theorem tells us how to choose between the two hypotheses to maximize P D for a given P FA . Based on the measurement model specified by Equation (13), the PDFs for the hypothesis H i (i = 0, 1) is given by: where Owing to the lack of knowledge about the signal s k , we cannot completely specify the PDFs under the hypothesis H 0 and apply the Likelihood Ratio Test (LRT) derived from Neyman-Pearson theorem. However, by replacing the unknown signal s k with its Maximum Likelihood Estimates (MLE), we can derive a Generalized Likelihood Ratio Test (GLRT). The MLE is obtained by maximizing Equation (16) with respect to the unknown parameters. Under the hypothesis H 0 , the MLEŝ k of the unknown s k can be obtained depending on the theory of MLE for vector parameters as follows: Then we can get:ŝ Substitutingŝ k into Equation (16) yields The GLRT decides on H 1 , if where the threshold γ is determined from If we combine Equations (18), (21) and (22), T G can be described as Simplifying yields Taking the natural logarithm on both sides of Equation (22) and simplifying yields: where γ = −2 ln γ. If the inequality in Equation (26) is satisfied, then the GLRT chooses the hypothesis that the magnetic field is quasi-static.

Heading Error Estimation
If the magnetic field is detected as quasi-static, the magnetically derived heading will be adopted to estimate the heading errors of the gyroscope based on the ZUPT-aided EKF algorithm. Let ω t be three-axis gyroscope output at time t, C t−1 denotes the rotation matrix that transforms from the body to the navigation frame at time t − 1, which is updated with the gyroscopic information. A Padé approximation of the exponential function is utilized to update orientation [18]: where ∆t denotes the sample interval, I 3 represents the unit matrix of three order, Ω t is the skew symmetric matrix of angular rates: Let δx t|t be the error state vector of EKF at time t: The vector δx t|t contains 5 components, namely the estimated errors in attitude (δφ t ), positions (δr t ), velocity (δv t ), the bias of the gyroscope (δω t ) and the accelerometer (δa t ). All these 5 components have three elements each, corresponding to the three-dimension estimation. The linearized state transition model is: where δx t|t−1 is the one step predicted error state, δx t−1|t−1 is the last filtered error state, w t−1 is the process noise with the covariance matrix Q t = E(w t w T t ). Φ t is the 15 × 15 state transition matrix: where 0 3 denotes 3 × 3 zero matrix and S(a t ) is the skew symmetric matrix for accelerations that allows the EKF to act as an inclinometer S(a t ) = a t = [a x t , a y t , a z t ] is the acceleration that has been transformed to the navigation frame. The measurement model is where z t is the error measurements, H is the measurement matrix, n t is the measurement noise with the covariance matrix G t = E(n t n T t ). If the magnetic field is detected as quasi-static during the stance phase, the error measurement m t can be calculated as: where v t is the velocity error derived from the accelerometer. ∆θ t is the heading error derived from the gyroscope: Here, as the magnetic field is detected as quasi-static, the heading θ m t derived from the magnetic field can be utilized as a reliable measurement to eliminate heading errors of the gyroscope. The measurement matrix H, for ZUPT update, is given by: Then the filtered error state δx t|t is obtained: where K t is the Kalman gain matrix. The attitude refinement is achieved by updating the rotation matrix C t with the attitude errors δφ t . Assuming δφ t is small, the corrected rotation matrix C t is given by: where Θ t is the skew symmetric matrix for small angles: According to the above description, we summarize the procedure of the proposed method and the flowchart is shown in Figure 4. According to the above description, we summarize the procedure of the proposed method and the flowchart is shown in Figure 4.

Results and Discussion
To analyze the performance of the proposed method, experiments have been carried out. As shown in Figure 5, the foot-mounted module we use in these experiments is a Multiple Inertial Measurement Units (MIMU) platform. The platform holds 8 MPU9250 IMUs (InvenSense Inc., San

Results and Discussion
To analyze the performance of the proposed method, experiments have been carried out. As shown in Figure 5, the foot-mounted module we use in these experiments is a Multiple Inertial Measurement Units (MIMU) platform. The platform holds 8 MPU9250 IMUs (InvenSense Inc., San Jose, CA, USA), with 4 on the top side and 4 on the bottom side. The data is stored and can be transferred to a computer via a USB interface. Each MPU9250 IMU consists of a three-axis MEMS accelerometer, a three-axis MEMS gyroscope and a three-axis MEMS magnetometer. The advantage of the MIMU platform has already been elaborated in Reference [32], so it will not be discussed here. The system is mounted in the heel of shoes and the sampling rate is 400 Hz.

Results and Discussion
To analyze the performance of the proposed method, experiments have been carried out. As shown in Figure 5, the foot-mounted module we use in these experiments is a Multiple Inertial Measurement Units (MIMU) platform. The platform holds 8 MPU9250 IMUs (InvenSense Inc., San Jose, CA, USA), with 4 on the top side and 4 on the bottom side. The data is stored and can be transferred to a computer via a USB interface. Each MPU9250 IMU consists of a three-axis MEMS accelerometer, a three-axis MEMS gyroscope and a three-axis MEMS magnetometer. The advantage of the MIMU platform has already been elaborated in Reference [32], so it will not be discussed here. The system is mounted in the heel of shoes and the sampling rate is 400 Hz.

Magnetometer Calibration Experiments
We conducted one test to assess the performance of the proposed magnetometer calibration method. In the test, the proposed method was compared with the CLSEF method. We performed the CLSEF method by rotating the sensor in all possible directions in a constant uniform magnetic field. The sensor measurement of rotation is depicted in Figure 6 and the calibration parameters of the two methods are summarized in Table 1. The variable c K denotes the combination of the soft iron errors, misalignment and scale factor errors and c B represent hard iron errors.

Magnetometer Calibration Experiments
We conducted one test to assess the performance of the proposed magnetometer calibration method. In the test, the proposed method was compared with the CLSEF method. We performed the CLSEF method by rotating the sensor in all possible directions in a constant uniform magnetic field. The sensor measurement of rotation is depicted in Figure 6 and the calibration parameters of the two methods are summarized in Table 1. The variable K c denotes the combination of the soft iron errors, misalignment and scale factor errors and B c represent hard iron errors.    The results for heading calibration by the two methods are compared in Figure 7. As the drift of gyroscope is very slow, it can be ignored during the magnetometer calibration procedure. Therefore, the heading derived from gyroscopes can be utilized as the reference to assess the performance of the two methods. In Figure 7, the difference between the reference and the proposed is no more than 5 • and much smaller than the CLSEF method. Although the CLSEF calibration method can provide accuracy calibration parameters in free space scenario, its performance will be degraded in the foot-mounted application because of the effects caused by the floor and shoes. The proposed calibration method can eliminate most of the effects and performs better in the foot-mounted application scenario.
The results for heading calibration by the two methods are compared in Figure 7. As the drift of gyroscope is very slow, it can be ignored during the magnetometer calibration procedure. Therefore, the heading derived from gyroscopes can be utilized as the reference to assess the performance of the two methods. In Figure 7, the difference between the reference and the proposed is no more than 5° and much smaller than the CLSEF method. Although the CLSEF calibration method can provide accuracy calibration parameters in free space scenario, its performance will be degraded in the footmounted application because of the effects caused by the floor and shoes. The proposed calibration method can eliminate most of the effects and performs better in the foot-mounted application scenario.

Heading Error Correction Experiments
To use the proposed detector effectively, there are several tuning parameters, denoted in Equation (26), need to be quantified.σ θ , σ I and σ B give a measure of the noise which is encountered during quasi-static field periods. To evaluate these parameters, the magnetometer and gyroscope data has been collected for more than two hours. The corresponding three kinds of noise variance are measured as σ θ = 0.067 rad/s, σ I = 0.054 rad/s and σ B = 5.35 mGs. The threshold γ can be selected based on the relationship between the probability of detection and the acceptable probability of false alarms. The Receiver Operating Characteristics (ROC) curve is utilized for this and it also defines the performance of the detector. The ROC of the detectors was calculated by varying the thresholds of the detectors and comparing the decisions of the detectors to the output of the reference system. We calculated the test statistics for the detectors using a window size N = 3. Figure 8 shows the ROC curve of the proposed detector. When test threshold is set to 170, the corresponding probability of detection is about 0.8. When the probability of detection is larger than 0.8, the false alarm probability increases much faster than the detection probability. Thus, we choose 170 as the test threshold.
probability of false alarms. The Receiver Operating Characteristics (ROC) curve is utilized for this and it also defines the performance of the detector. The ROC of the detectors was calculated by varying the thresholds of the detectors and comparing the decisions of the detectors to the output of the reference system. We calculated the test statistics for the detectors using a window size = 3 N . Figure 8 shows the ROC curve of the proposed detector. When test threshold is set to 170, the corresponding probability of detection is about 0.8. When the probability of detection is larger than 0.8, the false alarm probability increases much faster than the detection probability. Thus, we choose 170 as the test threshold. To evaluate the positioning performance of the proposed method, we conducted two tests in real indoor environments which are very common for pedestrian navigation application: office and shopping mall. In the office test, we walked along the corridors and walls in a typical office at normal speed (approximately 5 km/h) and then kept walking three loops with the same path. The total distance of the path was approximately 1018 m in 786 steps. Figure 9 gives the bird's eye view of the test region selected in the school where the authors work. To evaluate the positioning performance of the proposed method, we conducted two tests in real indoor environments which are very common for pedestrian navigation application: office and shopping mall. In the office test, we walked along the corridors and walls in a typical office at normal speed (approximately 5 km/h) and then kept walking three loops with the same path. The total distance of the path was approximately 1018 m in 786 steps. Figure 9 gives the bird's eye view of the test region selected in the school where the authors work.
probability of false alarms. The Receiver Operating Characteristics (ROC) curve is utilized for this and it also defines the performance of the detector. The ROC of the detectors was calculated by varying the thresholds of the detectors and comparing the decisions of the detectors to the output of the reference system. We calculated the test statistics for the detectors using a window size = 3 N . Figure 8 shows the ROC curve of the proposed detector. When test threshold is set to 170, the corresponding probability of detection is about 0.8. When the probability of detection is larger than 0.8, the false alarm probability increases much faster than the detection probability. Thus, we choose 170 as the test threshold. To evaluate the positioning performance of the proposed method, we conducted two tests in real indoor environments which are very common for pedestrian navigation application: office and shopping mall. In the office test, we walked along the corridors and walls in a typical office at normal speed (approximately 5 km/h) and then kept walking three loops with the same path. The total distance of the path was approximately 1018 m in 786 steps. Figure 9 gives the bird's eye view of the test region selected in the school where the authors work.   Figure 10 depicts the step-wise trajectories obtained using the pure ZUPT, iQSF methods and the proposed method, respectively. In Figure 10a the heading error accumulates due to gyro biases which cannot be observed using the pure ZUPT. The final positioning errors and heading errors (difference between the initial and final heading) reach up to 2 m and 14 • respectively. The iQSF method and the proposed method perform much better than the pure ZUPT algorithm as shown in Figure 10b,c. Both methods achieve positioning error of no more than 0.5 m. As most of the paths are straight lines in this test, so the iQSF method can also detect most potential quasi-static magnetic field and then estimates the heading errors of the gyroscope. Figure 10 depicts the step-wise trajectories obtained using the pure ZUPT, iQSF methods and the proposed method, respectively. In Figure 10a the heading error accumulates due to gyro biases which cannot be observed using the pure ZUPT. The final positioning errors and heading errors (difference between the initial and final heading) reach up to 2 m and 14° respectively. The iQSF method and the proposed method perform much better than the pure ZUPT algorithm as shown in Figure 10b,c. Both methods achieve positioning error of no more than 0.5 m. As most of the paths are straight lines in this test, so the iQSF method can also detect most potential quasi-static magnetic field and then estimates the heading errors of the gyroscope. Step-wise trajectory of the iQSF method; (c) Step-wise trajectory of the proposed method. Figure 11 shows the heading results of the first loop using the proposed method. It can be seen that although the magnetically derived heading sometimes changes violently due to the perturbations, there is the moment when it is stable and similar to the heading derived from gyroscope, which can be used to estimate the heading error originating from the gyroscope drift. Step-wise trajectory of the iQSF method; (c) Step-wise trajectory of the proposed method. Figure 11 shows the heading results of the first loop using the proposed method. It can be seen that although the magnetically derived heading sometimes changes violently due to the perturbations, there is the moment when it is stable and similar to the heading derived from gyroscope, which can be used to estimate the heading error originating from the gyroscope drift.
(c) Figure 10. The step-wise trajectories of the office test: (a) Step-wise trajectory of the pure ZUPT method; (b) Step-wise trajectory of the iQSF method; (c) Step-wise trajectory of the proposed method. Figure 11 shows the heading results of the first loop using the proposed method. It can be seen that although the magnetically derived heading sometimes changes violently due to the perturbations, there is the moment when it is stable and similar to the heading derived from gyroscope, which can be used to estimate the heading error originating from the gyroscope drift. Figure 11. The heading results of the first loop using the proposed method.
We conducted the second trajectory test in a shopping mall to evaluate the effectiveness of the proposed method in large-scale scenario. The total distance of this trajectory was approximately 2480 m in 1920 steps and 14 landmarks were set along the path to calculate the positioning errors. Figure  12 shows the trajectories obtained using pure ZUPT, iQSF and the proposed method. The maximum trajectory error obtained using the iQSF method is approximately 51 m whereas that for the proposed method is approximately 3 m. We conducted the second trajectory test in a shopping mall to evaluate the effectiveness of the proposed method in large-scale scenario. The total distance of this trajectory was approximately 2480 m in 1920 steps and 14 landmarks were set along the path to calculate the positioning errors. Figure 12 shows the trajectories obtained using pure ZUPT, iQSF and the proposed method. The maximum trajectory error obtained using the iQSF method is approximately 51 m whereas that for the proposed method is approximately 3 m.  Figure 12. Trajectories obtained using pure ZUPT, iQSF and the proposed method. Figures 13 and 14 depict the detection results of quasi-static magnetic field using iQSF and the proposed method, respectively. Compared with the office test, the trajectory of the shopping mall test is non-straight for most of the time. In this case, most of the potential quasi-static magnetic field cannot be detected using the iQSF method, while the proposed method is not subject to the straight line constraint. The quasi-static magnetic field can be detected successfully using the proposed method and then the gyroscope biases are estimated and compensated during the quasi-static magnetic field periods.  Figures 13 and 14 depict the detection results of quasi-static magnetic field using iQSF and the proposed method, respectively. Compared with the office test, the trajectory of the shopping mall test is non-straight for most of the time. In this case, most of the potential quasi-static magnetic field cannot be detected using the iQSF method, while the proposed method is not subject to the straight line constraint. The quasi-static magnetic field can be detected successfully using the proposed method and then the gyroscope biases are estimated and compensated during the quasi-static magnetic field periods. Figure 12. Trajectories obtained using pure ZUPT, iQSF and the proposed method. Figures 13 and 14 depict the detection results of quasi-static magnetic field using iQSF and the proposed method, respectively. Compared with the office test, the trajectory of the shopping mall test is non-straight for most of the time. In this case, most of the potential quasi-static magnetic field cannot be detected using the iQSF method, while the proposed method is not subject to the straight line constraint. The quasi-static magnetic field can be detected successfully using the proposed method and then the gyroscope biases are estimated and compensated during the quasi-static magnetic field periods.   During the test, we crossed the 14 landmarks for a total of 42 times and counted the positioning errors of the iQSF method and the proposed method as shown in Figures 15 and 16. The mean error trajectory obtained using the iQSF method has a mean positioning error of 11.3 m as compared to 1.5 m in case of the proposed method. Thus, the positioning errors are reduced by 85% by utilizing the proposed method for constraining the heading drift. It is worth mentioning here that all the test results are concluded from the proposed method without any aiding maps or other external signals. During the test, we crossed the 14 landmarks for a total of 42 times and counted the positioning errors of the iQSF method and the proposed method as shown in Figures 15 and 16. The mean error trajectory obtained using the iQSF method has a mean positioning error of 11.3 m as compared to 1.5 m in case of the proposed method. Thus, the positioning errors are reduced by 85% by utilizing the proposed method for constraining the heading drift. It is worth mentioning here that all the test results are concluded from the proposed method without any aiding maps or other external signals. errors of the iQSF method and the proposed method as shown in Figures 15 and 16. The mean error trajectory obtained using the iQSF method has a mean positioning error of 11.3 m as compared to 1.5 m in case of the proposed method. Thus, the positioning errors are reduced by 85% by utilizing the proposed method for constraining the heading drift. It is worth mentioning here that all the test results are concluded from the proposed method without any aiding maps or other external signals.

Conclusion
In this paper, a heading error estimation method for using the quasi-static magnetic field has been developed. The proposed approach is implemented in two sections. Firstly, to eliminate the errors originating from the system platform, floor and shoes, a magnetometer calibration method is proposed. In the proposed calibration method, the three-axis magnetometer calibration is transformed to that of two-axis according to the motion model of the pedestrian, hence the calibration procedure is simplified. Secondly, a novel quasi-static magnetic field detection approach is developed and utilized to suppress the heading drift of the gyroscope. The proposed method detects the quasi-static magnetic field using the change rate of the magnetic field magnitude, the inclination angle and the difference between the magnetically derived heading and the heading derived from gyroscope. The detected quasi-static magnetic field is utilized as the measurement fed into the ZUPTaided EKF to estimate the heading errors.
Selection of two real indoor environments (office and shopping mall) provided data sets for a detailed analysis of the performance of the proposed method. The experimental results indicate that the proposed method is capable of effectively suppressing the heading errors caused by the drift of gyroscopes, reducing the overall positioning error budget by over 85% in shopping mall environment. No other sensors of any kind are used in the proposed method, except IMUs and magnetometers. The approach developed herein are self-contained and assumed a denied GNSS environment. However, GNSS is partly available in numerous indoor environments and urban

Conclusion
In this paper, a heading error estimation method for using the quasi-static magnetic field has been developed. The proposed approach is implemented in two sections. Firstly, to eliminate the errors originating from the system platform, floor and shoes, a magnetometer calibration method is proposed. In the proposed calibration method, the three-axis magnetometer calibration is transformed to that of two-axis according to the motion model of the pedestrian, hence the calibration procedure is simplified. Secondly, a novel quasi-static magnetic field detection approach is developed and utilized to suppress the heading drift of the gyroscope. The proposed method detects the quasi-static magnetic field using the change rate of the magnetic field magnitude, the inclination angle and the difference between the magnetically derived heading and the heading derived from gyroscope. The detected quasi-static magnetic field is utilized as the measurement fed into the ZUPT-aided EKF to estimate the heading errors.
Selection of two real indoor environments (office and shopping mall) provided data sets for a detailed analysis of the performance of the proposed method. The experimental results indicate that the proposed method is capable of effectively suppressing the heading errors caused by the drift of gyroscopes, reducing the overall positioning error budget by over 85% in shopping mall environment.
No other sensors of any kind are used in the proposed method, except IMUs and magnetometers. The approach developed herein are self-contained and assumed a denied GNSS environment. However, GNSS is partly available in numerous indoor environments and urban canyons. Therefore, we will integrate the two approaches to maximize availability and accuracy in the future.