An Aeromagnetic Compensation Method Based on a Multimodel for Mitigating Multicollinearity.

Aeromagnetic surveys play an important role in geophysical exploration and many other fields. In many applications, magnetometers are installed aboard an aircraft to survey large areas. Due to its composition, an aircraft has its own magnetic field, which degrades the reliability of the measurements, and thus a technique (named aeromagnetic compensation) that reduces the magnetic interference field effect is required. Commonly, based on the Tolles-Lawson model, this issue is solved as a linear regression problem. However, multicollinearity, which refers to the case when more than two model variables are highly linearly related, creates accuracy problems when estimating the model coefficients. The analysis in this study indicates that the variables that cause multicollinearity are related to the flight heading. To take this point into account, a multimodel compensation method is proposed. By selecting the variables that contribute less to the multicollinearity, different sub-models are built to describe the magnetic interference of the aircraft when flying in different orientations. This method restricts the impact of multicollinearity and improves the reliability of the measurements. Compared with the existing methods, the proposed method reduces the interference field more effectively, which is verified by a set of airborne tests.

1. Introduction

Aeromagnetic surveys originated in the 1930s for military applications [1] and now play a very important role in many other fields, such as geophysical exploration. With the magnetometers installed aboard aircrafts, people can measure the magnetic field at a very flexible scale [2]. However, since the aircraft has its own magnetic field, which degrades the reliability of the measurements, a technique to reduce the magnetic interference field is required, that is, aeromagnetic compensation [3,4].

Aeromagnetic compensation is mainly based on the Tolles–Lawson (T–L) model [5,6], which decomposes the magnetic interference field into three sources: the permanent field, induced field, and eddy-current field. Combining all of these fields, the magnetic interference field can be described as a linear equation with 18 or 16 terms (the latter being simplified from an 18-term equation) [7]. Solving this equation is called a calibration, which is a key issue in aeromagnetic compensation, and the elements of the solution are called coefficients. Commonly, a figure-of-merit (FOM) flight [8] is implemented for calibration, which includes four orthogonal headings and three sets of maneuvers (pitches, rolls, and yaws) in each heading. After removing the ambient magnetic field from the measured total field, the remaining field is considered to be the magnetic interference field, which can be the dependent variable of the 16-term (or 18-term) linear equation. Then the coefficients can be estimated through regression.

To yield a more accurate set of coefficients, researchers have made many efforts toward improving the model [9,10], optimizing the solving method [11,12,13,14,15,16], and correcting the sensor errors [17]. Theoretically, with the assumption that the ambient magnetic field is uniform [7], the coefficients of the T–L model are directly linked to the aircraft itself because they are due to its properties, such as the materials and the electrical systems. However, in practice, the coefficients are hard to obtain accurately. One reason for this is that the uniformity assumption of the ambient magnetic field is unrealistic. Another significant factor is multicollinearity among the 16 or 18 variables of the linear equation, which causes noise sensitivity in the estimation. To mitigate the multicollinearity, some statistical methods are utilized, such as ridge regression (RR) [7], truncated singular value decomposition (TSVD) [14], and partial least-squares regression (PLSR) [15]. These methods render the estimated coefficients more accurate for the 18-term model but are not always effective for the 16-term model [7]. Two variables are excluded in the 16-term model because they can be linearly represented by other variables and contribute significantly to the multicollinearity. Nevertheless, the multicollinearity is typically still strong.

In this paper addressing aeromagnetic compensation based on scalar magnetometers, we analyze the sources of multicollinearity and find those that depend on the flight heading. Differing from the present methods that regard the FOM flight as a whole, here we propose a multimodel method to compensate for the magnetic interference field of the aircraft, according to the flight heading. By selecting different variable sets for different headings, multicollinearity can be inhibited.

This paper is structured as follows. In Section 2, we describe the T–L model, analyze why multicollinearity occurs, and present our method. In Section 3, the method’s performance is verified by a set of airborne tests. Section 4 is the conclusion.

2. Analysis and Method

2.1. T–L Model

To better describe the analysis and the method, the T–L model is briefly introduced here. First, the reference system is displayed in Figure 1. The reference system is attached to the aircraft: the origin O is at the location of the magnetometer, the L-axis is parallel to the left, the T-axis points toward the nose, and the V-axis is vertical. The vector represents the magnetic field of the earth, and X, Y, and Z are the angles between the magnetic field of the earth and the three axes, respectively. The vector represents the total field measured by the magnetometer. As scalar magnetometers only measure the intensity of the total magnetic field and as the intensity of the earth magnetic field is a few orders greater than that of the interference field, we can define the scalar of the interference field as the projection from itself onto the direction of the earth’s magnetic field.

Figure 1

The reference system.

According to the T–L model [6], the magnetic interference field can be decomposed into three parts: the permanent field, the induced field, and the eddy-currents field. The permanent field is constant, written as , the scalar of which can be represented as

The second part of the magnetic interference field is induced by the magnetic field of the earth and can be represented as

The third part is caused by the eddy currents in the aircraft and can be expressed as where , , and are the time derivatives of , , and , respectively, which are defined as where W represents X, Y, or Z, which varies with time. Thus, the magnetic interference is expressed as

This is the 18-term T–L equation.

As we have the term can be represented by the and terms, and so it can be removed. Furthermore, considering the derivatives of the two sides of (6), namely the term can be removed too. Lastly, the 16-term T–L model is built.

2.2. Analysis of Multicollinearity

Section 2.1 indicates that the variables of the T–L model are formed by the direction cosines of the magnetic field of the earth. These cosines can be calculated from the aircraft attitude. Let and represent the magnetic heading and the dip angle, respectively, and , , and represent the pitch, roll, and yaw angles, respectively. Then, the direction cosines of the magnetic field of the earth can be represented as below [6]. For pitches:

For rolls:

For yaws:

As is a certain value in each heading, some simplified expressions can be deduced from Equations (8)–(10). In Table 1 and Table 2, we list the model variables and the relevant simplified expressions when the heading is north and east, respectively. The cases when the heading is south and west are very similar to Table 1 and Table 2, respectively.

Table 1

Model variables in the north heading.

Variable	Expression	Pitch	Roll	Yaw
A_1	cosX	0	sinϕsinψ	−cosϕsinΩ
A_2	cosY	cosϕcosλ+sinϕsinλ	cosϕ	cosϕcosΩ
A_3	cosZ	sinϕcosλ−cosϕsinλ	sinϕcosψ	sinϕ
A_4	cos^2X	0	sin^2ϕsin^2ψ	cos^2ϕsin^2Ω
A_5	cosXcosY	0	sinϕcosϕsinψ	−cos^2ϕsinΩcosΩ
A_6	cosXcosZ	0	sin^2ϕsinψcosψ	−sinϕcosϕsinΩ
A_7	cos^2Y	cos^2ϕcos^2λ+2sinϕcosϕsinλcosλ+sin^2ϕsin^2λ	cos^2ϕ	cos^2ϕcos^2Ω
A_8	cosYcosZ	sinϕcosϕcos^2λ+(sin^2ϕ−cos^2ϕ)sinλcosλ−sinϕcosϕsin^2λ	sinϕcosϕcosψ	sinϕcosϕcosΩ
A_9	cosXcosX˙	0	sin^2ϕsinψcosψ	cos^2ϕsinΩcosΩ
A_10	cosXcosY˙	0	0	cos^2ϕsin^2Ω
A_11	cosXcosZ˙	0	−sin^2ϕsin^2ψ	0
A_12	cosYcosX˙	0	sinϕcosϕcosψ	−cos^2ϕcos^2Ω
A_13	cosYcosY˙	sinϕcosϕcos^2λ+(sin^2ϕ−cos^2ϕ)sinλcosλ−sinϕcosϕsin^2λ	0	−cos^2ϕsinΩcosΩ
A_14	cosYcosZ˙	−cos^2ϕcos^2λ−2sinϕcosϕsinλcosλ−sin^2ϕsin^2λ	−sinϕcosϕsinψ	0
A_15	cosZcosX˙	0	sin^2ϕcos^2ψ	−sinϕcosϕcosΩ
A_16	cosZcosY˙	sin^2ϕcos^2λ−2sinϕcosϕsinλcosλ+cos^2ϕsin^2λ	0	−sinϕcosϕsinΩ

Table 2

Model variables in the east heading.

Variable	Expression	Pitch	Roll	Yaw
A_1	cosX	cosϕ	cosϕcosψ+sinϕsinψ	cosϕcosΩ
A_2	cosY	sinϕsinλ	0	cosϕsinΩ
A_3	cosZ	sinϕcosλ	sinϕcosψ−cosϕsinψ	sinϕ
A_4	cos^2X	cos^2ϕ	cos^2ϕcos^2ψ+2sinϕcosϕsinψcosψ+sin^2ϕsin^2ψ	cos^2ϕcos^2Ω
A_5	cosXcosY	sinϕcosϕsinλ	0	cos^2ϕsinΩcosΩ
A_6	cosXcosZ	sinϕcosϕcosλ	sinϕcosϕcos^2ψ+(sin^2ϕ−cos^2ϕ)sinψcosψ−sinϕcosϕsin^2ψ	sinϕcosϕcosΩ
A_7	cos^2Y	sin^2ϕsin^2λ	0	cos^2ϕsin^2Ω
A_8	cosYcosZ	sin^2ϕsinλcosλ	0	sinϕcosϕsinΩ
A_9	cosXcosX˙	0	sinϕcosϕcos^2ψ+(sin^2ϕ−cos^2ϕ)sinψcosψ−sinϕcosϕsin^2ψ	−cos^2ϕsinΩcosΩ
A_10	cosXcosY˙	sinϕcosϕcosλ	0	cos^2ϕcos^2Ω
A_11	cosXcosZ˙	−sinϕcosϕsinλ	−cos^2ϕcos^2ψ−2sinϕcosϕsinψcosψ−sin^2ϕsin^2ψ	0
A_12	cosYcosX˙	0	0	−cos^2ϕsin^2Ω
A_13	cosYcosY˙	sin^2ϕsinλcosλ	0	cos^2ϕsinΩcosΩ
A_14	cosYcosZ˙	−sin^2ϕsin^2λ	0	0
A_15	cosZcosX˙	0	sin^2ϕcos^2ψ−2sinϕcosϕsinψcosψ+cos^2ϕsin^2ψ	−sinϕcosϕsinΩ
A_16	cosZcosY˙	sin^2ϕcos^2λ	0	sinϕcosϕcosΩ

As , , and are small angles, the following assumptions are established: and

According to Equations (11) and (12), we have some linear correlations in Table 1:

These equations indicate that when the heading is the north, , , , and correlate with each other, and , , and also correlate with each other. These correlations lead to multicollinearity. For the T–L model, the correlations mean that some coefficients can be represented by other coefficients. Hence, these variables are ineffective in calculating the magnetic interference field and negative for compensation.

When the heading is east, multicollinearity also exists; however, the correlations are caused by different sources. From Table 2, the linear correlations are as below: where , , , and correlate with each other, and , , and correlate with each other.

The comparison of these two cases indicates that the variables causing the multicollinearity are different. In other words, the variables causing multicollinearity are related to the flight heading.

2.3. The Multimodel Method

Based on the theory of regression, multicollinearity renders the variances of the least squares (LS) estimates large, and setting some coefficients to zero can improve the estimate accuracy [18]. Thus, reducing some of the correlated variables can be effective in mitigating the multicollinearity of the T–L model. However, in Section 2.2, we found that the correlated variables are not the same in different headings. So treating the different heading cases in different ways can suppress the influence of the multicollinearity more effectively. In this paper, we propose a method to select a subset of the variables to build a sub-model for a specific heading.

We choose the variance inflation factors (VIFs), which are indicators of multicollinearity [19] and defined as the diagonal elements of , to be the criteria for determining the selected variables. For the purpose of retaining the variables that contribute less to multicollinearity, the selection manner is designed to exclude the variables that have larger VIFs because the larger the VIF is, the more the relevant variable contributes to the multicollinearity and should not be selected. If the variables that have larger VIFs are excluded, the multicollinearity of the shrunk model will be weakened.

The implementation procedures of our method are shown in Figure 2 and described as follows. In the calibration stage, we divide the raw data collected during a calibration flight into four segments, according to the flight headings. For each segment, we select a variable set, by comparing the VIFs, to build the sub-model and estimate the coefficients. The estimation method is still the LS method. From the calibration flight, four sets of coefficients are obtained. In the compensation stage, first, the current heading needs to be identified to choose the relevant set of coefficients, and then the compensation is implemented.

Figure 2

Procedures for calibration and compensation. VIF—variance inflation factor.

Note that in the calibration flight, large maneuvers cause extra displacements of the magnetometer. Due to the gradient of the magnetic field of the earth, the measurement of the magnetic field will be also disturbed by these displacements, besides the aircraft’s magnetic field. To abandon the impact caused by the displacement of the magnetometer in the magnetic field of the earth, the gradient should be corrected. According to [7], we add a three-term gradient model after each shrunk T–L model. If the earth’s magnetic gradient is written as , where , , and represent the gradient components of the magnetic north, magnetic east, and the local vertical, respectively, the variation in the magnetic field caused by the displacement of the magnetometer can be represented by where , , and represent the magnetometer displacements in the magnetic north, the magnetic east, and the local vertical, respectively. Finally, the sub-model corresponding to the specific flight heading is the combination of the shrunk T–L model and the gradient model.

3. Experiment

3.1. Implementation

The experiment was implemented in a Y-12 aircraft with the magnetometers installed in the tail boom. As shown in Figure 3, the scalar magnetometer and the three-axis vector magnetometer were equipped at the pole and the middle of the tail boom, respectively. The scalar magnetometer was a homemade helium magnetometer with sensitivity ; its principle can be seen in [20,21].

Figure 3

The aircraft used in the experiment.

As shown on the right of Figure 3, two datasets were adopted to test the compensation performance of the proposed method. They were labeled Dataset-A and Dataset-B. The flight headings were sorted by north, east, south, and west. The maneuvers were operated in the order of roll, pitch, and yaw. The flight altitude was about 3000 m. In each heading, the flight length was about 2 to 3 min.

Table 3 shows the process of selecting the variables in the north heading case, where FVS is short for full variable set and EVS is short for excluded variable set. In the beginning, we entered all of the 16 variables (FVS) and listed their VIFs. The largest VIF belonged to , namely . In other words, contributed the most to the multicollinearity and should be excluded. Subsequently, the other variables (EVS-1) were entered. This time, many of the VIFs were much smaller because the relevant variables exhibited strong collinearity with . The largest VIF belonged to (), which was also excluded. In EVS-2, the excluded variable was () and in EVS-3 it was (). After this, all of the VIFs were much smaller.

Table 3

Variance inflation factors (VIFs) in the north heading. FVS—full variable set; EVS—excluded variable set.

Variables			VIFs
	FVS	EVS-1	EVS-2	EVS-3	EVS-4
A_1	9730.40	6640.77	5119.77	—	—
A_2	159,672.01	—	—	—	—
A_3	112,862.35	39,651.18	23.55	23.54	23.32
A_4	16,679.02	9875.18	21.07	20.81	20.78
A_5	1907.04	1312.74	1311.42	512.57	511.42
A_6	4727.15	3772.57	2369.45	509.71	507.73
A_7	36,281.72	35,521.73	—	—	—
A_8	26,334.11	526.00	20.67	20.67	20.49
A_9	626.94	626.33	623.44	619.19	15.07
A_10	1.60	1.60	1.60	1.46	1.43
A_11	2.51	2.49	2.47	2.46	2.35
A_12	488.07	488.07	487.06	482.17	443.24
A_13	1946.64	1945.84	1932.61	1913.33	—
A_14	1169.43	1168.18	1162.95	1154.38	25.10
A_15	483.08	483.08	481.86	477.28	436.85
A_16	510.18	510.11	504.93	499.07	10.90

The east case was similar to the north case. Table 4 shows the VIFs of the east segment in the same dataset. In the beginning, although the VIF of () was larger, we retained and excluded the second largest, (); this was because and are already excluded in the 16-term model, and so deleting causes the loss of too much information about the Z-axis. By deleting and , the VIF of is small. Subsequently, (), (), and () were excluded, identically to the excluded variables in the north case but with the former Y replaced by X. The X-axis points north if the aircraft is oriented to the east, instead of the Y axis. These results indicate that the excluded variable set is stable.

Table 4

Variance inflation factors (VIFs) in the east heading.

Variables			VIFs
	FVS	EVS-1	EVS-2	EVS-3	EVS-4
A_1	397,320.73	—	—	—	—
A_2	27,811.26	27,266.38	–	—	—
A_3	(400,873.05)	2809.42	2763.78	3.38	3.38
A_4	35,405.72	3384.67	3336.12	—	—
A_5	3626.27	3534.16	37.41	37.39	37.29
A_6	51,546.73	49.68	49.51	3.37	3.36
A_7	138.65	3.30	3.03	2.61	1.79
A_8	11,488.87	11,299.04	39.48	39.41	39.09
A_9	328.74	324.86	222.25	219.56	—
A_10	71.27	67.20	66.60	66.60	66.16
A_11	104.81	104.04	75.54	73.93	4.41
A_12	22.37	22.35	19.98	18.54	18.47
A_13	1.60	1.60	1.49	1.45	1.39
A_14	19.68	19.67	16.97	15.74	15.72
A_15	98.89	97.12	61.96	61.87	5.20
A_16	69.11	65.23	64.58	64.36	63.81

In some data segments, the case where causes the largest VIF happens sometimes. However, to avoid losing too much information, we retain and exclude the variable that causes the second-largest VIF. Taking no account of this case, the results of Table 3 and Table 4 are stable. In the south and west segments, the subsets are the same as in the north and east segments, respectively. In other datasets, these results are reiterated. With these subsets, we build a specific model for each orientation.

Figure 4 shows the compensation results of our proposed method, where the selected subset was EVS-4; the coefficients calculated from Dataset-A were used to compensate Dataset-B. The FOMs, which are defined as the sums of peak-to-peak values in 12 maneuvers, were 26.6769 nT and 7.0414 nT before and after compensation, respectively. The standard deviations (STDs) of the uncompensated and compensated signals were 0.5230 nT and 0.1187 nT. The improvement ratio (IR), which is defined as the ratio between these two STDs, was 4.4054.

Figure 4

Compensation results.

3.2. Comparison with Conventional Methods

To verify the performance of the multimodel (MM) method proposed in this paper, we compare it with two conventional methods: the traditional LS method and the RR method. The three methods were respectively tested on the two datasets. Figure 5 shows the comparison results, where the FOMs and IRs of the three methods are exhibited. In this figure, the labels of the x-axis indicate the relationships between the two datasets. The first letter represents the calibration dataset, and the second one represents the compensation dataset (e.g., AB means using the coefficient set calculated from Dataset-A to compensate Dataset-B). From Figure 5, we find that in each test, the LS method and the RR method had the worst performances, and the MM method showed an obvious improvement.

Figure 5

Comparisons between the three compensation methods: (a) figure-of-merit (FOM), (b) standard deviation (STD), and (c) improvement ratio (IR). LS—least squares; RR—ridge regression; MM—multimodel.

Figure 6 shows the compensation results of the three methods in the four flight headings of Dataset-B, where the coefficient sets were calculated from Dataset-A. The figure indicates that the conventional LS and RR methods, based on the 16-term model, gave similar results, and the proposed MM method performed better in compensating for the magnetic interference field. The peak-to-peak values are showed in Table 5. In each box, the four numbers represent the peak-to-peak values of uncompensated, LS-compensated, RR-compensated, and MM-compensated signals, respectively. For these methods, The improvement factors of the FOMs were 2.1192, 2.1448, and 3.7886, respectively.

Figure 6

Compensation results of the three compensation methods. In the legends, represents the standard deviation (STD), with unit of nT.

Table 5

The peak-to-peak values of uncompensated (UN), LS-compensated, RR-compensated, and MM-compensated signals (UN|LS|RR|MM).

Maneuver	North	East	South	West	Sum
Roll	1.19|1.39|1.40|0.72	5.11|1.33|1.41|0.61	2.01|1.29|1.26|0.34	3.88|2.00|1.89|0.86	12.19|6.00|5.97|2.53
Pitch	2.71|0.66|0.62|0.47	0.89|0.34|0.32|0.36	3.04|0.24|0.23|0.15	1.83|1.86|1.85|1.74	8.47|3.09|3.02|2.73
Yaw	1.25|0.67|0.65|0.31	1.84|1.03|1.06|0.42	1.10|0.69|0.67|0.17	1.82|1.11|1.07|0.88	6.01|3.49|3.45|1.78
Sum	5.15|2.72|2.67|1.50	7.84|2.70|2.80|1.39	6.15|2.21|2.16|0.66	7.54|4.96|4.81|3.79	26.68|12.59|12.44|7.04

3.3. Compensations in Level Flights

The FOM flights mentioned above, containing different headings and different large maneuvers, are implemented for calibrating the coefficients and evaluating the compensation ability. Here, we present two tests on level flights, which are closer to the practical conditions. These flights, named Line-A and Line-B, were oriented to the North and South, respectively, with a length of 3 min. The coefficients applied here were those calculated from Dataset-A. Figure 7 shows the compensation results. The STD of Line-A reduced from 0.0506 nT to 0.0206 nT and that of Line-B reduced from 0.1076 nT to 0.0452 nT. The IRs were 2.4565 and 2.3835, respectively. This figure illustrates that the proposed MM method can help to measure the magnetic field more accurately with an airborne magnetometer.

Figure 7

Compensation results: (a) Line-A and (b) Line-B.

We also compared the MM method with the two conventional methods on Line-A and Line-B. Table 6 shows the results of the comparisons between the three methods. From this, we find that on both datasets, the proposed MM method gave better results than the other two conventional methods.

Table 6

The STDs and IRs of the three methods on Line-A and Line-B.

Dataset		STD (nT)			IR
	LS	RR	MM	LS	RR	MM
Line-A	0.0223	0.0216	0.0206	2.2695	2.3452	2.4565
Line-B	0.0579	0.0594	0.0452	1.8594	1.8124	2.3835

3.4. Compensation in Non-Standard Headings

The experimental results above have proved that in the standard headings (north, east, south, and west), the MM method is effective in aeromagnetic compensation and better than the conventional methods based on the 16-term T–L model. However, in practice, sometimes the flight heading is not standard. Hence, it is necessary to test the performance of the proposed MM method in non-standard headings.

Figure 8 shows the tracks of other datasets, labeled as Dataset-C and Dataset-D. These two datasets were collected while the aircraft was flown towards the Northeast, Southeast, Southwest, and Northwest. The aircraft, the magnetometer, and the other sensors were unchanged. The compensation result is shown in Figure 9. The FOM was reduced from 28.3546 nT (uncompensated) to 4.4953 nT (compensated), the STD was reduced from 0.5813 nT (uncompensated) to 0.0902 nT (compensated), and the IR was 6.4470.

Figure 8

The tracks of Dataset-C and Dataset-D.

Figure 9

The compensated signal of Dataset-C. The coefficients are calculated from Dataset-C.

We also compared our MM method with the LS and RR methods. The comparison results are shown in Figure 10. It indicates that the MM method is still better than the traditional LS and RR methods in the non-standard flight headings.

Figure 10

Comparisons between the three compensation methods in non-standard headings: (a) figure-of-merit (FOM), (b) standard deviation (STD), and (c) improvement ratio (IR).

4. Conclusions

Magnetometers are usually equipped aboard an aircraft in many applications, such as geophysical exploration. As the magnetic field caused by the aircraft interferes with the measurements, an aeromagnetic compensation method should be applied. Most aeromagnetic compensation methods are based on the T–L model and are restricted by its multicollinearity. Herein, we have found that the model variables causing the multicollinearity differ according to the flight heading, based on which we proposed a multimodel method to mitigate the multicollinearity. This method built different sub-models for different headings by selecting the variables with smaller VIFs. In the real flight experiments, the MM method reduced the FOM from 26.6769 to 7.0414. The improvement factor is 3.7886, higher than the factors yielded by two conventional methods (the LS and RR methods), which were 2.1192 and 2.1448, respectively. In the level flight tests, the MM method reduced the STDs by about 2.4 times.