An Optimized Two-Step Magnetic Correction Strategy by Means of a Lagrange Multiplier Estimator with an Ellipsoid Constraint.

The geomagnetic field is as fundamental a constituent of passive navigation as Earth's gravity. In cases where no other external attitude reference is available, for the direct heading angle estimation by a typical magnetic compass, a two-step optimized correction algorithm is proposed to correct the model coefficients caused by hard and soft iron nearby. Specifically, in Step 1, a Levenberg-Marquardt (L-M) fitting estimator with an ellipsoid constraint is applied to solve the hard magnetic coefficients. In Step 2, a Lagrange multiplier estimator is used to deal with the soft magnetic iron circumstance. The essential attribute of "the two-step" lies in its eliminating the coupling effects of hard and soft magnetic fields, and their mutual interferences on the pure geomagnetic field. Under the conditions of non-deterministic magnetic interference sources with noise, the numerical simulation by referring to International Geomagnetic Reference Field (IGRF), and the laboratory tests based upon the turntable experiments with Honeywell HMR3000 compass (Honeywell, Morristown, NJ, USA) conducted, the experimental results indicate that, in the presence of the variation of multi-magnetic interferences, the RMSE (Root Mean Square Error) value of the estimated total magnetic flux density by the proposed two-step estimator falls to 0.125 μT from its initial 2.503 μT, and the mean values of the heading angle error estimates are less than 1°. The proposed solution therefore, exhibits ideal convergent properties, fairly meeting the accuracy requirements of non-tactical level navigation applications.

1. Introduction

Magnetic compasses are usually considered effective measuring units for heading angle estimates, the accuracy of which, in general, will determine the performances of the whole navigation system [1,2]. In practical measurements, the measured data derived from uncontrolled environmental magnetic fields will lead to large error in evaluating the relevant attitude parameters [3,4]. Among the error sources, hard iron (mainly refers to the permanent magnets, the magnetized iron or steel) and soft iron (mainly refers to the iron which is easily magnetized once being placed in a magnetic field environment, especially the pure iron) interferences around the carriers, are generally the main error sources for magnetic measuring results [5]. Therefore, before being put into use, it is essential to correct the orthogonal tri-axis measurements of magnetic compasses. The research on dealing with magnetic correction can be classified into three categories. One is to utilize the reference information derived from, for example, Global Positioning System (GPS), the inertial measurement units (IMUs), etc., to correct the direct measuring results of magnetic compasses [6]. One is to use the external devices to measure the orthogonal tri-axis attitude angles of magnetic compasses as assisted sources [7]. Another one is to fulfill the magnetic correction by a class of algorithm-based means, termed ellipsoid fitting solutions [8].

Note that, differing from the hardware-upgrading schemes stated in the first two categories, for the ellipsoid fitting solutions, there is no need to take into account the cost, sizes or power consumption of the external devices. Therefore, for practical applications it is relatively feasible to adopt the algorithm-based means that concentrates on algebraic ellipsoid fitting and magnetic coefficient estimates. Since the direct solutions involve complex coupling consequences, the preferred solution therefore, is intended to underpin progressive estimation in many domains such as simplifying the calculating iteration, eliminating the mutual coupling effects, and reducing the computational work [9,10,11]. As a typical instance, D. Gebre-egziabher [12] proposed a nonlinear estimation algorithm to calibrate the solid-state strapdown magnetometers, first putting forward a two-step model. In Step 1, they skillfully introduced a group of auxiliary variables to denote nine parameters in the typical ellipsoid constraint, linearizing an equation set for the nine parameter estimates rather than performing the direct calculation for the hard and soft magnetic coefficients. In Step 2, these nine parameters were solved algebraically. Experiments showed that this two step estimator outperformed an extended Kalman filter formulation [12,13].

Inspired by this, in this framework, a proposed two-step magnetic correction algorithm with fusing the L-M fitting concept and the Lagrange multiplier idea is expected to eliminate the negative influences the hard and soft iron cause. The essence of which consists in the algebraic coefficient estimates corresponding to an irregular ellipsoid to a normal sphere transform by two steps. In view of the analyses above, the optimized correction algorithm appears fairly reliable, being characterized by better independence of the coefficient estimates in the presence of mutual coupling effects between hard and soft magnetic field interferences.

The outline of the remainder of the paper is as follows: in the following section, the preprocessing of the magnetic data (more precisely, signals) is first stated. In Section 3, the optimized correction for hard magnetic coefficients by the L-M fitting estimator is described, and the Lagrange multiplier method process is also detailed allowing for ideal precision of the soft magnetic coefficient estimates. The numerical simulation and laboratory tests are carried out in Section 4, which emphasizes the experimental assessment on geomagnetic field measurement error estimates extracted from a calibrated compass under the conditions of various non-deterministic interferences in terms of, location, sample point shape, form (hard or soft iron) and intensity. Section 5 presents the main conclusion of this investigation.

2. Preprocessing of Magnetic Signals

Affected by internal structural and external environmental magnetic field uncertainties, the magnetic compass measurements are subject to noise interferences. In view of the frequency domain, improving the signal-to-noise ratio (SNR) of raw signals appears to be a principal solution to high-precision navigation service on the basis of external sensors’ measurements [14]. In these frameworks, our work mainly consists of three parts:

We adopted the empirical mode decomposition (EMD) method to reduce the noise and filter the raw signals [15,16], and we obtained a series of multiple intrinsic mode function (IMF) components of the measured signals;

For the noisy IMF components, we adopted the improved wavelet threshold denoising method (WTD) for signal noise reduction, and this process can be seen as completely retaining the pure IMF components;

We further reconstructed a series of desired signals by fusing the original pure IMF components and the de-noised IMF components.

Since our experiments involve numerical simulation and actual laboratory tests, we emphasize that, for the generated data (the simulated magnetic flux density data), the Gaussian white noise is successfully eliminated by the above EMD theory. Analogously, in the real-time laboratory measuring situations, the smoothing requirements are also fully satisfied.

3. Two-Step Optimized Magnetic Coefficients Correction Algorithm

When tilting a magnetic compass in any direction, observation shows that the shape of trajectory derived from orthogonal tri-axis outputs with respect to carrier coordinate reference maintains a normal sphere, which is characterized by the hypothesis that the magnetic compass is working in a pure geomagnetic environment without any interference like hard or soft iron on the host platform [17,18]. In practice, however, even though the intrinsic instrumentation errors of the magnetic compass, for example, scale factors, non-orthogonality, and offsets, exert certain influences on measuring results, we concentrate more upon eliminating the impacts of environmental magnetic field interferences, because in cases where hard or soft iron exists, we can only obtain an irregular ellipsoid rather than an approximate normal sphere [19,20,21]. The following subsections describe the elaborate tri-axis distribution transformation by mathematical means.

3.1. 1st Step Correction for Hard Magnetic Coefficients with L-M Fitting Estimator

Let the carrier coordinate reference be built along three axes of the magnetic compass. The hard and soft iron interference matrixes are, respectively, denoted as H and [22]. With real and measuring value matrix of geomagnetic field being denoted by H and H respectively, the error compensation model of compass is obtained as follows [23]:

Let , rearranging Equation (1):

Let , then Equation (2) can be further written as:

Since matrix M satisfies the Jo Riski decomposition [24], so M is symmetric and it follows that: where a, b,…, f denote the unknown elements of matrix M. On the assumption of , and , Equation (3), therefore, satisfies the following expression:

According to Equation (5), for n sets of H, we get n linear equations in terms of the form below: where , .

Let (i = 1, 2, …, 9) indicate the elements of known vector k (being evaluated beforehand by L-M fitting method), the hard iron interference matrix with respect to k, therefore, can be calculated from the ternary equations set below:

3.2. 2nd Step Correction for Soft Magnetic Coefficients with Lagrange Multiplier Estimator

Let being a collection of data points after hard iron interference compensation, so that Equation (3) reduces to:

Equation (8) in this case can be rewritten by substituting Equation (4), it follows that:

Equation (9) represents an ellipsoid with the following constraint: where, , , q is a unknown positive integer.

For each , define a new vector in augmented form, , where i = 1, 2, …, N representing the measuring set number.

To estimate the coefficients a, b, c, d, e and f in Equation (9), let be a coefficient vector in general augmented form, whose values are to be determined.

Thus, rewrite Equation (9) in the following matrix form, and clearly it can be expanded to a set of linear equations: where D = [ X’ X’ … X’’], . Denote the geometric distance (the minimum distance between the measured and the optimal fitted ellipsoid) by the matrix below:

Minimizing Equation (12) under the constraint indicated in Equation (10) is our direct objective. To achieve this, rewrite Equation (10) in terms of , and we get:

The robust solution to the above optimization problem, therefore, turns to a matter for skilled mathematic iteration.

Construct the simultaneous equations by typical Lagrange multiplier model: where is the scalar Lagrange multiplier, and represents “Lagrange function of”.

Differentiate the Lagrange function with respective of and respectively, and let and , we obtain:

Note here, Equation (15) has only one solution when q > 3, which denotes the general eigenvector associated with the unique positive eigenvalue of a generalized eigenvalue equation .

Solving for Equations (15) and (16), we reconstruct and in the following forms: where, the matrixes P11, P12, P22 are of size 6 × 6, 6 × 1, and 1 × 1 respectively. Clearly, ,.

Hence we are dealing with an eigensystem of two equations given by: where, extracting the elements located in the first 6th row and 6th column of original matrix . In Equation (20), P22 would be nonsingular in case of the relevant sampling points being not coplanar (since groups of data are derived from a 3-dimensional ellipsoid), then it can be verified that Equation (20) is solvable, which yields:

Substituting Equation (21) into Equation (19), we obtain the generalized eigenvalue system equation [25], it follows that:

As stated above, is no less than 3. When , , then the quadrique indicated by Equation (5) must be an ellipsoid, we say that is the sufficient condition for the present ellipsoid. In this case, we need to determine a bigger integer value , so as to obtain a wider range of ellipsoid quadrique and satisfy .

A rule of thumb, satisfactory in simple cases, is to designate big enough (like 50), then we locate by means of bisecting, viz., iteratively calculating it in the intervals between and .

The algorithm flow chart for the geomagnetic measurement error correction is shown in Figure 1. As indicated, the implementing volume has been divided into two parts, and the coefficients estimation appears to be merely forward, by means of which, we say the coupling effects’ elimination experienced by hard and soft magnetic field interferences is fairly achievable. For detailed experiments and result analyses, please refer to Section 4.

Figure 1

Algorithm flow chart for the geomagnetic measurement error correction.

4. Experiments

Experimental observations consist of numerical simulations and laboratory tests. We are generally considering the simulation environment as the ideal geomagnetic field [26], since the pure geomagnetic field never exists under normal conditions of sequent tests.

4.1. Numerical Simualtions and Analyses

The geomagnetic field of the Jilin area (Jilin City, China) is designated as our simulation environment and the geomagnetic field distribution information is derived from the International Geomagnetic Reference Field (IGRF 12) by reference, with altitude 183.4 m in the urban area, latitude 43°52′ N, longitude 126°33′ E, and date June 2018 [27]. Two sets of values for the soft magnetic coefficient matrix and hard magnetic coefficient matrix H are designated, as denoted by HH (benchmark value) in Table 1. Similarly, the standard deviation of the white noise (μT) is set to be 0.1 and 0.01, respectively.

Table 1

Simulation results.

HH	ε_3×3	[0.83574−0.046350.0265601.125640.01152001.10233]		[0.935520.096540.0523401.146210.03207001.19657]
	Hi	[0.90.7−0.4]		[0.3-0.30.8]
σ(μT)		0.1	0.01	0.1	0.01
q		7	8	9	9
HH1	[a]Hi	[0.90070.6986−0.3990]	[0.90030.6988−0.4010]	[0.2991-0.30080.7983]	[0.2992-0.29980.8006]
	[b]Hi	[0.90350.6162−0.3780]	[0.90440.6861−0.3941]	[0.3035-0.30620.9092]	[0.2150-0.31010.7390]
	[a]ε_3×3	[0.83601−0.046020.0267501.125370.01180001.10221]	[0.83594−0.046110.0267801.125410.01175001.10229]	[0.935960.096120.0526501.146020.03176001.19640]	[0.936010.096040.0525401.145920.03169001.19636]
	[b]ε_3×3	[0.83635−0.047110.0269601.125930.01211001.10286]	[0.83629−0.047050.0269501.125890.01208001.10281]	[0.936250.097870.0515601.145610.03286001.19689]	[0.936060.097210.0516101.145600.03279001.19687]
HP(μT)	[a]	0.847	0.074	0.831	0.092
	[b]	0.941	0.085	0.970	0.183
HM(°)	[a]	0.656	0.230	0.696	0.247
	[b]	0.764	0.315	0.765	0.317
HT(s)	[a]	1.05	1.02	1.09	0.99
	[b]	1.42	1.35	1.44	1.39

The simulated magnetic flux density in μTesla is derived from a Matlab simulation platform, and the preprocessing work focuses on smoothing the white noise involved. On the basis of the de-noising data, the orthogonal tri-axis magnetic field distribution before and after correction is illustrated by the dots in Figure 2, it is actually easier to indicate that the shape of dots changes from an irregular ellipsoid to an approximate normal sphere.

Figure 2

Tri-axis magnetic field distribution.

Figure 3 shows the error of total magnetic flux density. As presented, the error with raw data is very large as predicted (see Figure 3a), and the result (0.7 μT) by the proposed two-step estimator (see Figure 3c) exhibits better than that (1 μT) by the traditional two-step estimator stated in Reference [12,13] (see Figure 3b). It should be noted that the average of the total magnetic flux density conforms with the pre-set condition reference, i.e., B=54.397 μT (see Figure 2), and the solution adopted has been proved to possess strong convergence.

Figure 3

Total magnetic flux density before and after correction. (a) the error with raw data; (b) the result by the traditional two-step estimator; (c) the result by the proposed two-step estimator.

More detailed simulation results are given in Table 1, the derived q is based on a typical bisecting iteration. HH1 represents the estimated hard and soft magnetic coefficients, here [a] and [b] respectively represent the results by the proposed two-step estimator and the traditional two-step estimator, HP (μT) represents the maximum error of total magnetic flux density, HM (°) represents the mean of the heading angle error and HT (s) represents the average running time for the numerical simulation under the same conditions. The simulation results show that the soft and hard magnetic coefficients derived by the L-M fitting and the Lagrange multiplier estimator are much closer to the true values. In contrast, the optimized algorithm-based scheme shortened the executing time, leading to better heading angle estimation accuracy without implementing any hardware-upgrading means or integrating any external measurement units.

4.2. Laboratory Tests and Evaluations

The correction is experimentally assessed to evaluate the effectiveness of the optimized solution in practice. The laboratory tests are conducted based upon the turntable experiments with Honeywell HMR3000 compass (which measurement accuracy is within 100 nT) and the real geomagnetic field reference is provided by a high-precision proton magnetometer (the measurement accuracy is within 1 nT). To begin with, the compass is mounted on a stationary base on a two-axis horizontal turntable. In sequence, in the course of collecting orthogonal 3-axis magnetic flux density data, the compass rotates at a constant speed to capture uniformly distributed samples across different directions. Three experiments are implemented indicating different locations of interference sources, as illustrated in Figure 4. Specifically, it is to ensure that the hard and soft magnetic material be distributed around the magnetic compass in concentric circles. Referring to Reference [28], the radii of three interference sources are set to be 90, 70 and 50 cm, respectively.

Figure 4

Interference source distribution.

As in the numerical simulation, the preprocessing of raw magnetic data should first be examined. Taking Y-axis measurements as an example, Figure 5a presents raw magnetic data derived from the Honeywell HMR3000 compass. Since the raw magnetic data contains noise, in order to filter the noise, the raw data was decomposed into six IMF components by the EMD method, diagrammatically represented by Figure 5b. In sequence, denoising the first three IMF components by improved WTD and fusing them with the remaining three IMF components, we achieved the final de-noised magnetic data. Figure 5c indicates the de-noised magnetic data by signal reconstruction in the frequency domain. Clearly, the de-noised data appears smoother, being better suited and favorable for the following ellipsoid fitting and relevant magnetic coefficients estimation.

Figure 5

The preprocessing of raw measuring magnetic data.

Figure 6 shows the fitted ellipsoids in cases of different circumstances where the orthogonal tri-axis magnetic field intensity is almost evenly distributed, evenly distributed, and partially distributed. Obviously, even the geomagnetic field distribution is more likely to form a regular ellipsoid, leading to an approximate normal sphere. In practice, however, it may well be that the turntable fails to ideally rotate at a constant speed, and the measured data consequently covers a very small portion of the ellipsoid surface.

Figure 6

The orthogonal tri-axis magnetic field distribution and the fitted ellipsoid.

In this case, determining the accurate radian and radial length at the very start is of great importance, which is achieved through multiple iterations of q. In our framework, the proposed two-step optimized algorithm by the L-M fitting and the Lagrange multiplier estimator has been applied to a vastly wider range of datasets, and the fitting results are adequately satisfied.

The estimated total geomagnetic field measurement error is diagrammatically represented in Figure 7, with a-c corresponding to different interference sources nearby, referring to Figure 4, the radii of the interference sources appear to be 90 cm, 70 cm and 50 cm respectively. Clearly, the traditional two-step estimator reveals its limited capacity in more accurate error estimation, and the undesired curve wave also illustrates degrees of instability. The algorithm stability still needs to be enhanced. By comparison, the proposed two-step estimator presents the optimal error estimation irrespective of the variation of the environmental magnetic fields, whose results exhibit better convergent properties with respect to the coupling effects experienced by hard and soft iron.

Figure 7

The estimated total magnetic flux density error. (a) the interference source radius is 90 cm; (b) the interference source radius is 70 cm; (c) the interference source radius is 50 cm.

The corresponding stochastic data (root mean square error, RMSE) comparison is illustrated by the bar graph in Figure 8. As shown, for the case that the hard and soft iron are very close to the magnetic compasses, the RMSE value of the estimated total magnetic flux density extracted from the proposed two-step estimator falls to 0.125 μT from the initial 2.503 μT (with raw data from Honeywell HMR3000). Specifically, the traditional two-step estimator presents a relatively large RMSE, and this is due to the fact that in step 2, in order to calculate the relevant nine parameters algebraically, a number of multivariate equations are involved with hardly guaranteeing the independence of the equations, that is, the solutions for nine parameters would be not unique, and thus the magnetic flux density error would be large. In contrast, even though the proposed two-step estimator follows the same line as the traditional one in the two-step magnetic correction design, the optimized estimator addresses this problem by means of multiple iterations of the equations involved, and after these proactive iterations, the calculation error decreases, leading to more ideal RMSE results.

Figure 8

The comparison of RMSE. (a) the interference source radius is 90 cm; (b) the interference source radius is 70 cm; (c) the interference source radius is 50 cm.

The estimated heading angle error before and after correction is presented in Figure 9. Clearly, whose peak values even approach 20° when directly using the raw data from Honeywell HMR3000, i.e., the measured data is not corrected at all (see Figure 9a). It should be noted that, in this case, the distance of the interference object nearby is set to be 90 cm, indicating that the interferences are not too severe. Through Matlab simulative calculations, the mean values of the estimated heading angle error corresponding to Figure 9b–d (whose radii of interference sources are 90, 70 and 50 cm, respectively) are 0.486°, 0.825° and 0.918°, respectively. It is obvious to indicate that, after the magnetic correction, the mean values of the heading angle error estimates are no more than 1° regardless of the difference of locations of interference sources and radii of concentric circles, and that the optimized two-step magnetic correction strategy fairly meets the accuracy requirements of non-tactical level navigation applications.

Figure 9

The estimated heading angle error before and after correction by the proposed two-step estimator. (a) before correction; (b) the interference source radius is 90 cm; (c) the interference source radius is 70 cm; (d) the interference source radius is 50 cm.

5. Conclusions

An optimized two-step model by the L-M fitting and the Lagrange multiplier estimator for magnetic field measurement correction was employed, after exhibiting ideal convergent properties in the numerical simulations, which was experimentally assessed with the Honeywell HMR3000 compass to eliminate the magnetic interferences caused by surrounding hard and soft iron. It demonstrated that, compared with the traditional two-step estimator, the proposed estimator effectively reduces the environmental magnetic field disturbance and its accuracy and reliability remain optimal with the variation of non-deterministic interferences. The experimental results indicate that the mean values of the heading angle error estimates are less than 1° under strong magnetic interferences. The algorithm-based solution with no extra assistant measurement units, therefore, may be considered as a reference tool for a class of magnetic field coefficients correction with environmental interferences hypothesis.