An Improved Aeromagnetic Compensation Method Robust to Geomagnetic Gradient

: Aeromagnetic surveys play an important role in many ﬁelds, for example, archaeology, anti-submarine warfare, and geophysical exploration. Being in the geomagnetic ﬁeld, the aircraft generates a great deal of magnetic interference, resulting in bad performance during detection surveys. Thus, it is necessary and important to perform aeromagnetic compensation in advance. Conventional aeromagnetic compensation methods consider that the geomagnetic gradient is approximately zero after bandpass ﬁltering, bringing about the inaccuracy of compensation coefﬁcients. To address this issue, an improved aeromagnetic compensation method robust to geomagnetic gradient is proposed. In this study, the International Geomagnetic Reference Field (IGRF) model was employed to model the geomagnetic gradient. Then, the estimated geomagnetic gradient was subtracted from the measured data, which improved the accuracy of the compensation equations. Field experiments were conducted to verify the effectiveness of the proposed method. The experimental results show that compared to the traditional method, the compensation performance of the proposed method was improved by 152% to 329%. For the level ﬂight, the standard deviation of residual noise after compensation can be as low as 3.3pT. The results indicate that the proposed method can signiﬁcantly improve the compensation effect, showing great beneﬁts for weak magnetic anomaly detection.


Introduction
Aeromagnetic surveying has been widely used in geophysical exploration, unexploded ordnance (UXO) detection, anti-submarine warfare, and many other fields [1][2][3][4][5][6][7]. The magnetic anomaly signal is measured by the magnetometer installed on the flight platform. However, the detection performance of the aeromagnetic survey system is significantly affected by the magnetic interference field that is generated by the flight platform. Therefore, it is important to estimate and eliminate the interference field [8][9][10][11], which method is called compensation.
The earliest aeromagnetic compensation model was established by Tolles [12,13]. The magnetic interference was divided into three categories, namely the permanent, the induced, and the eddy current field, respectively. In 1961, Leliak [14] proposed an effective calibration flight scheme and derived a method for solving the compensation coefficients. Since then, Leach [15] applied the least square method to solve the compensation coefficients. The multicollinearity in the compensation model was analyzed, and the ridge experimental results indicated that the accuracy of the magnetometer in the yaw angle estimation had been greatly enhanced compared to the traditional ellipsoid fitting method.
The total magnetic field, measured by an optically pumped magnetometer (OPM), includes the geomagnetic field and the magnetic interference field. For the traditional compensation method, the geomagnetic gradient is usually considered a constant value, which is negligent after bandpass filtering. However, the geomagnetic field is a function of the spatial position and time, as described by the International Geomagnetic Reference Field (IGRF) model [27]. Based on the latest 13th generation IGRF model published by the International Association of Geomagnetism and Aeronomy (IAGA), the global geomagnetic intensity is shown in Figure 1. According to the requirements of standard calibration flight, the distance of profiles in each direction ranges from a few kilometers to more than ten kilometers, resulting in a geomagnetic gradient of several nanoTeslas or even hundreds of nanoTeslas. For the high-precision magnetic survey system, this is interference that cannot be ignored. Some previous studies have proved that considering the influence of the geomagnetic gradient can improve the accuracy of estimated compensation coefficients. Dou et al. [19] applied the World Magnetic Model 2010 to a generate geomagnetic field, which was then used to solve the compensation coefficients. The simulation results showed that the compensation performance can be improved. Ge et al. [28] proposed a relatively simple geomagnetic model to solve the compensation coefficients. The flight test results showed that the compensation performance of the two flights was improved by 27% and 22%, respectively. Zhou et al. [29] analyzed the influence of the geomagnetic gradient on the towed bird system and established the towed bird interference model; the feasibility of this compensation method was verified by actual flight tests. To eliminate the residual interference field caused by the geomagnetic gradient as far as possible, an improved aeromagnetic compensation method robust to geomagnetic gradient is proposed in this paper. The application steps are as follows: first, the high-precision flight profiles are obtained by the Differential Global Positioning System (DGPS); secondly, the geomagnetic field information is obtained according to the IGRF model; finally, the background field, including the geomagnetic gradient, is substituted into the compensation model to obtain the compensation coefficients. Flight tests were carried out, and the results show that the residual noise of the proposed method was lower than that of the traditional method, which demonstrates that the geomagnetic gradient was effectively compensated; meanwhile, a higher detection ability was obtained for the aeromagnetic survey system.
The detailed study results are presented in the following four sections: Section 2 summarizes the principle of the traditional Tolles-Lawson model and the calibration flight. In Section 3, a new aeromagnetic compensation model robust to the geomagnetic gradient is extended by introducing the IGRF model, which was subsequently used to solve the compensation coefficients. An aeromagnetic survey system with low noise instruments and its data analysis are presented in Section 4. Section 5 presents the conclusion of this paper.

Tolles-Lawson Model
To compensate for the platform-generated interference field, a three-axis fluxgate magnetometer is used to measure the attitude of the platform according to the classical Tolles-Lawson (T-L) model [12,13]. The reference coordinate system attached to the aircraft is established, as shown in Figure 2, where the origin O represents the location of the magnetometer, and the T, L, and V are the transverse, longitudinal, and vertical axis, respectively. He represents the geomagnetic field, and X, Y, and Z are the angles between the geomagnetic field and the three axes, respectively. To eliminate the residual interference field caused by the geomagnetic gradient as far as possible, an improved aeromagnetic compensation method robust to geomagnetic gradient is proposed in this paper. The application steps are as follows: first, the high-precision flight profiles are obtained by the Differential Global Positioning System (DGPS); secondly, the geomagnetic field information is obtained according to the IGRF model; finally, the background field, including the geomagnetic gradient, is substituted into the compensation model to obtain the compensation coefficients. Flight tests were carried out, and the results show that the residual noise of the proposed method was lower than that of the traditional method, which demonstrates that the geomagnetic gradient was effectively compensated; meanwhile, a higher detection ability was obtained for the aeromagnetic survey system.
The detailed study results are presented in the following four sections: Section 2 summarizes the principle of the traditional Tolles-Lawson model and the calibration flight. In Section 3, a new aeromagnetic compensation model robust to the geomagnetic gradient is extended by introducing the IGRF model, which was subsequently used to solve the compensation coefficients. An aeromagnetic survey system with low noise instruments and its data analysis are presented in Section 4. Section 5 presents the conclusion of this paper.

Tolles-Lawson Model
To compensate for the platform-generated interference field, a three-axis fluxgate magnetometer is used to measure the attitude of the platform according to the classical Tolles-Lawson (T-L) model [12,13]. The reference coordinate system attached to the aircraft is established, as shown in Figure 2, where the origin O represents the location of the magnetometer, and the T, L, and V are the transverse, longitudinal, and vertical axis, respectively. H e represents the geomagnetic field, and X, Y, and Z are the angles between the geomagnetic field and the three axes, respectively. The direction cosines of the angles between the geomagnetic field and the platform are denoted as 1 , 2 , and 3 , and can be expressed as follows [21]: According to the T-L model, the platform-generated interference field can be grouped into three types, namely the permanent, the induced, and the eddy current magnetic fields. Define Hd(t) as the total intensity of the interference field, then Hd(t) can be expressed as follows: The permanent magnetic field HPERM(t) caused by the permanent magnetism of the aircraft's ferromagnetic parts can be expressed as follows: where ci (i = 1, 2, 3) is the compensation coefficient of the permanent magnetic field, and Ai(t) (i = 1, 2, 3) is the variable constituted by the direction cosine. The induced magnetic field HIND(t) is created in paramagnetic parts by the geomagnetic field, and it can be expressed as follows: where ci (i = 4, …, 9) is the compensation coefficient of the induced magnetic field, and Ai(t) (i = 4, …, 9) is the variable constituted by the direction cosine and He(t).
The eddy current, generated by the maneuvers of the metal components of the aircraft, is proportional to the time change rate of the flux through these components. The eddy current interference field HEDDY(t) can be expressed as follows:  The direction cosines of the angles between the geomagnetic field and the platform are denoted as u 1 , u 2 , and u 3 , and can be expressed as follows [21]: According to the T-L model, the platform-generated interference field can be grouped into three types, namely the permanent, the induced, and the eddy current magnetic fields. Define H d (t) as the total intensity of the interference field, then H d (t) can be expressed as follows: The permanent magnetic field H PERM (t) caused by the permanent magnetism of the aircraft's ferromagnetic parts can be expressed as follows: where c i (i = 1, 2, 3) is the compensation coefficient of the permanent magnetic field, and A i (t) (i = 1, 2, 3) is the variable constituted by the direction cosine. The induced magnetic field H IND (t) is created in paramagnetic parts by the geomagnetic field, and it can be expressed as follows: where c i (i = 4, . . . , 9) is the compensation coefficient of the induced magnetic field, and A i (t) (i = 4, . . . , 9) is the variable constituted by the direction cosine and H e (t). The eddy current, generated by the maneuvers of the metal components of the aircraft, is proportional to the time change rate of the flux through these components. The eddy current interference field H EDDY (t) can be expressed as follows: where u 1 , u 2 , and u 3 are the time derivatives of u 1, u 2 , and u 3 , respectively, c i (i = 10, . . . , 18) is the compensation coefficient of the eddy current magnetic field, and A i (t) (i = 10, . . . , 18) is the variable constituted by the direction cosine and H e (t). The interference field generated by the platform can be expressed as follows: The corresponding matrix notation of Equation (6) can be expressed as follows: where H d and C are both column vectors constituted by H d (t) and c i (i = 1, 2, . . . , 18), and A is the feature matrix that can be expressed as follows: The compensation coefficients can be obtained by the least square algorithm [15,30], as given by: Under the Tolles-Lawson model, the compensation coefficients are only related to the platform and the local background geomagnetic field. Thus, the compensation coefficients do not change unless the airplane platform system has changed. The time-varying A i (t) (i = 1, 2, . . . , 18) can be obtained from the direction cosines u i (i = 1, 2, 3), their time derivatives u i (i = 1, 2, 3), and the magnitude of the geomagnetic field. It can be clearly demonstrated that the interference field generated by the platform would be well eliminated if the accurate compensation coefficients can be obtained.

Calibration Flight
To obtain the compensation parameters of the platform, the aircraft needs to conduct a set of maneuvers according to certain requirements, also known as calibration flight [14,31]. As shown in Figure 3, the calibration flight is carried out in four orthogonal headings, and the pitch, the yaw, and the roll maneuvers are completed successively in each direction, for which the standard amplitudes of the pitch, the roll, and the yaw are ±5 • , ±5 • , and ±10 • , respectively. The altitude of the calibration flight is usually set at several kilometers to avoid being affected by the magnetic interference from the ground. Usually, the calibration flight is followed by a test flight that is used to evaluate the effectiveness of the estimated compensation parameters. where u1′, u2′, and u3′ are the time derivatives of u1, u2, and u3, respectively, ci (i = 10, …, 18) is the compensation coefficient of the eddy current magnetic field, and Ai(t) (i = 10, …, 18) is the variable constituted by the direction cosine and He(t). The interference field generated by the platform can be expressed as follows: The corresponding matrix notation of Equation (6) can be expressed as follows: where Hd and C are both column vectors constituted by Hd(t) and ci (i = 1, 2, …, 18), and A is the feature matrix that can be expressed as follows: The compensation coefficients can be obtained by the least square algorithm [15,30], as given by: Under the Tolles-Lawson model, the compensation coefficients are only related to the platform and the local background geomagnetic field. Thus, the compensation coefficients do not change unless the airplane platform system has changed. The time-varying Ai(t) (i = 1, 2, …, 18) can be obtained from the direction cosines ui (i = 1, 2, 3), their time derivatives ui′ (i = 1, 2, 3), and the magnitude of the geomagnetic field. It can be clearly demonstrated that the interference field generated by the platform would be well eliminated if the accurate compensation coefficients can be obtained.

Calibration Flight
To obtain the compensation parameters of the platform, the aircraft needs to conduct a set of maneuvers according to certain requirements, also known as calibration flight [14,31]. As shown in Figure 3, the calibration flight is carried out in four orthogonal headings, and the pitch, the yaw, and the roll maneuvers are completed successively in each direction, for which the standard amplitudes of the pitch, the roll, and the yaw are ±5°, ±5°, and ±10°, respectively. The altitude of the calibration flight is usually set at several kilometers to avoid being affected by the magnetic interference from the ground. Usually, the calibration flight is followed by a test flight that is used to evaluate the effectiveness of the estimated compensation parameters.

Proposed Aeromagnetic Compensation Model
The total field intensity H T measured by the optically pumped magnetometer is the superposition of the geomagnetic field H e and the platform interference field H d , hence the platform interference field H d can be expressed as follows: It is difficult to obtain an accurate platform interference field in practice. A widely used method is to apply the bandpass filter to Equation (10), which can be expressed as follows: where bpf () is a bandpass finite impulse response (FIR) filter operator. The passband of the filter is determined by the frequency characteristics of the aircraft's maneuvers, which are measured by the three-axis fluxgate. The power spectrum of the direction cosines u 1, u 2 , and u 3 are estimated, as shown in Figure 4. It can be seen that the maneuvering frequency of the aircraft during the calibration flight is concentrated between 0.15 Hz and 0.5 Hz, which provides a reference for the design of the bandpass filter in data preprocessing. On the other hand, the frequency band of magnetic anomaly detection (MAD) is about 0.04 Hz-0.4 Hz. Therefore, the cut-off frequency of the bandpass filter is set between 0.04 Hz and 0.5 Hz.

Proposed Aeromagnetic Compensation Model
The total field intensity HT measured by the optically pumped magnetometer is the superposition of the geomagnetic field He and the platform interference field Hd, hence the platform interference field Hd can be expressed as follows: It is difficult to obtain an accurate platform interference field in practice. A widely used method is to apply the bandpass filter to Equation (10), which can be expressed as follows: where bpf() is a bandpass finite impulse response (FIR) filter operator. The passband of the filter is determined by the frequency characteristics of the aircraft's maneuvers, which are measured by the three-axis fluxgate. The power spectrum of the direction cosines u1, u2, and u3 are estimated, as shown in Figure 4. It can be seen that the maneuvering frequency of the aircraft during the calibration flight is concentrated between 0.15 Hz and 0.5 Hz, which provides a reference for the design of the bandpass filter in data preprocessing. On the other hand, the frequency band of magnetic anomaly detection (MAD) is about 0.04 Hz-0.4 Hz. Therefore, the cut-off frequency of the bandpass filter is set between 0.04 Hz and 0.5 Hz. For the traditional compensation method, the geomagnetic gradient is usually considered to be negligent, and the following equation can be obtained: However, for a high-precision aeromagnetic survey, it is critical to consider the influence of the geomagnetic gradient when estimating the compensation coefficients. As suggested by the IGRF model, the geomagnetic field can be expressed by the gradient of a scalar potential, and the potential function can be described as a finite series expansion of spherical harmonic coefficients, as follows: where r is the radial distance from the center of the Earth, θ and ϕ are the geocentric colatitude and longitude, respectively, t represents time, α = 6371.2 km is the approximate For the traditional compensation method, the geomagnetic gradient is usually considered to be negligent, and the following equation can be obtained: However, for a high-precision aeromagnetic survey, it is critical to consider the influence of the geomagnetic gradient when estimating the compensation coefficients. As suggested by the IGRF model, the geomagnetic field can be expressed by the gradient of a scalar potential, and the potential function can be described as a finite series expansion of spherical harmonic coefficients, as follows: where r is the radial distance from the center of the Earth, θ and φ are the geocentric co-latitude and longitude, respectively, t represents time, α = 6371.2 km is the approximate mean radius of the Earth. g n m (t) and h n m (t) are the Gauss coefficients, P n m (cos θ) is Schmidt semi-normalized associated Legendre functions of degree n and order m. N is the number of expansion terms of spherical harmonic function. The spherical coordinate system is defined in Figure 5. By deriving the scalar potential of the internal magnetic field, the zonal component X, meridional component Y, and radial component Z of the geomagnetic field in the spherical coordinate system can be obtained as follows: Appl. Sci. 2022, 12, 1490 7 of 16 mean radius of the Earth. gn m (t) and hn m (t) are the Gauss coefficients, Pn m (cos θ) is Schmidt semi-normalized associated Legendre functions of degree n and order m. N is the number of expansion terms of spherical harmonic function. The spherical coordinate system is defined in Figure 5. By deriving the scalar potential of the internal magnetic field, the zonal component X, meridional component Y, and radial component Z of the geomagnetic field in the spherical coordinate system can be obtained as follows: The geomagnetic field He can be expressed as follows: According to the aforementioned analysis, Equation (15) can be rewritten as a function related to the year, longitude, latitude, and height, as follows: where ψ is a model that can calculate the geomagnetic field by using the IGRF coefficients. Based on the magnetic field data from the magnetic survey satellites, observatories, and other magnetic survey programmers around the world, the IGRF coefficients are computed every five years by the IAGA. The 13th Generation IGRF, the latest version released in December 2019, represents the geomagnetic field from 2020 to 2025. The longitude, latitude, and height can be converted to the parameters r, ϕ, and θ in Equation (13) by coordinate transformation. Based on Equations (9), (11), and (16), a compensation method robust to geomagnetic gradient is proposed, and the following equation can be established: where HT = (HT(1), HT(2), …, HT(n)) T is the column vector consisting of the intensities of the total field, He = (He(1), He(2), …, He(n)) T is the column vector consisting of the geomagnetic field intensities, and bps(ℵ) represents the bandpass filter that operates on each column of the matrix ℵ. The classic least squares method can be then used to estimate the compensation coefficients. The geomagnetic field H e can be expressed as follows: According to the aforementioned analysis, Equation (15) can be rewritten as a function related to the year, longitude, latitude, and height, as follows: where ψ is a model that can calculate the geomagnetic field by using the IGRF coefficients. Based on the magnetic field data from the magnetic survey satellites, observatories, and other magnetic survey programmers around the world, the IGRF coefficients are computed every five years by the IAGA. The 13th Generation IGRF, the latest version released in December 2019, represents the geomagnetic field from 2020 to 2025. The longitude, latitude, and height can be converted to the parameters r, φ, and θ in Equation (13) by coordinate transformation. Based on Equations (9), (11), and (16), a compensation method robust to geomagnetic gradient is proposed, and the following equation can be established: where H T = (H T (1), H T (2), . . . , H T (n)) T is the column vector consisting of the intensities of the total field, H e = (H e (1), H e (2), . . . , H e (n)) T is the column vector consisting of the geomagnetic field intensities, and bps(ℵ) represents the bandpass filter that operates on each column of the matrix ℵ. The classic least squares method can be then used to estimate the compensation coefficients.

Aeromagnetic Survey System
As shown in Figure 6, the Cessna-208 fixed-wing aircraft was selected as the survey platform, and a 5-m-long probe rod was installed at the tail of the aircraft. The probe rod was equipped with an optically pumped cesium magnetometer and a three-axis fluxgate magnetometer. The data recording and compensation system was installed in the aircraft cabin. All these instruments were developed by the Aerospace Information Research Institute of the Chinese Academy of Sciences (AIRCAS). A DGPS antenna was attached to the top of the plane to better receive satellite signals.

Aeromagnetic Survey System
As shown in Figure 6, the Cessna-208 fixed-wing aircraft was selected as the survey platform, and a 5-m-long probe rod was installed at the tail of the aircraft. The probe rod was equipped with an optically pumped cesium magnetometer and a three-axis fluxgate magnetometer. The data recording and compensation system was installed in the aircraft cabin. All these instruments were developed by the Aerospace Information Research Institute of the Chinese Academy of Sciences (AIRCAS). A DGPS antenna was attached to the top of the plane to better receive satellite signals. The optically pumped cesium magnetometer was the most important sensor in the system, and the sensitivity was as low as 0.3 pT/√Hz at 1 Hz. The sensitivity of the fluxgate magnetometer was 6 pT/√Hz at 1 Hz. The technical specifications are shown in Tables 1  and 2.    The aeromagnetic data recording and compensation system and the Differential Global Positioning System (DGPS) receiver were installed in the cabin in Figure 7. A software application was used to check the status of sensors and DGPS. All of the devices were powered by lithium batteries. The technical specifications of the recording system and the DGPS are shown in Tables 3 and 4.     The block diagram of the aeromagnetic data recording and compensation system is shown in Figure 8. Optically pumped magnetometer data and three-axis magnetometer data were recorded on the rising edge of the pulse per second (PPS), which was obtained from the DGPS. The sensitivity of the frequency counter channel and the analog to digital converter channel were 0.1 pT/ √ Hz at 1 Hz and 1 pT/ √ Hz at 1 Hz, respectively. In the meantime, the longitude, latitude, and altitude were collected for aeromagnetic data compensation in real time. The raw data and compensated data were all recorded in the flight.
Appl. Sci. 2022, 12, 1490 10 of 16 The block diagram of the aeromagnetic data recording and compensation system is shown in Figure 8. Optically pumped magnetometer data and three-axis magnetometer data were recorded on the rising edge of the pulse per second (PPS), which was obtained from the DGPS. The sensitivity of the frequency counter channel and the analog to digital converter channel were 0.1 pT/√Hz at 1 Hz and 1 pT/√Hz at 1 Hz, respectively. In the meantime, the longitude, latitude, and altitude were collected for aeromagnetic data compensation in real time. The raw data and compensated data were all recorded in the flight.

Data Evaluation Methods
The key data evaluation methods are the standard deviations of compensated and uncompensated signals and the improvement ratio (IR) derived from them. The IR can be expressed as follows [32]: where σu and σc are, respectively, the standard deviation of the uncompensated and compensated signals. The larger the IR, the better the compensation result.

Compensation Flight and Level Flight
The compensation flight for this experiment was carried out in the South China Sea at an altitude of more than 3000 m and a flight speed of 180 km/h. Five cycles of maneuvers consisting of about ±5° rolls, ±5° pitches, and ±10° yaws are carried out in each heading. As it is shown in Figure 9, a calibration flight and a test flight were recorded, and these two flight paths were very similar. Usually, the test flight is conducted under the same magnetic background as the calibration flight to verify the performance of aeromagnetic compensation.

Data Evaluation Methods
The key data evaluation methods are the standard deviations of compensated and uncompensated signals and the improvement ratio (IR) derived from them. The IR can be expressed as follows [32]: where σ u and σ c are, respectively, the standard deviation of the uncompensated and compensated signals. The larger the IR, the better the compensation result.

Compensation Flight and Level Flight
The compensation flight for this experiment was carried out in the South China Sea at an altitude of more than 3000 m and a flight speed of 180 km/h. Five cycles of maneuvers consisting of about ±5 • rolls, ±5 • pitches, and ±10 • yaws are carried out in each heading. As it is shown in Figure 9, a calibration flight and a test flight were recorded, and these two flight paths were very similar. Usually, the test flight is conducted under the same magnetic background as the calibration flight to verify the performance of aeromagnetic compensation.
The compensation flight for this experiment was carried out in the South China Sea at an altitude of more than 3000 m and a flight speed of 180 km/h. Five cycles of maneuvers consisting of about ±5° rolls, ±5° pitches, and ±10° yaws are carried out in each heading. As it is shown in Figure 9, a calibration flight and a test flight were recorded, and these two flight paths were very similar. Usually, the test flight is conducted under the same magnetic background as the calibration flight to verify the performance of aeromagnetic compensation.  The level flight path is shown in Figure 10. The level flight was about 30 km, the aim of which was to verify the best performance of aeromagnetic compensation.
Appl. Sci. 2022, 12, 1490 11 of 16 The level flight path is shown in Figure 10. The level flight was about 30 km, the aim of which was to verify the best performance of aeromagnetic compensation.

Data Analysis
The compensation parameters were solved according to the process shown in Figure  11. The geomagnetic gradient was considered in the process of solving the compensation coefficients. The raw data were sampled at 160 Hz, then filtered by a lowpass filter with a cut-off frequency of 20 Hz. The geomagnetic field column He was computed with Equation (16), and the information matrix A was constituted by Equation (6). Finally, the compensation parameters were solved by Equation (9). After the calibration flight, the effectiveness of the proposed method was examined based on the test flight data.

Data Analysis
The compensation parameters were solved according to the process shown in Figure 11. The geomagnetic gradient was considered in the process of solving the compensation coefficients. The raw data were sampled at 160 Hz, then filtered by a lowpass filter with a cut-off frequency of 20 Hz. The geomagnetic field column H e was computed with Equation (16), and the information matrix A was constituted by Equation (6). Finally, the compensation parameters were solved by Equation (9). After the calibration flight, the effectiveness of the proposed method was examined based on the test flight data.
The compensation coefficients could be obtained without taking into consideration the geomagnetic gradient (i.e., the traditional method), as shown in Table 5. Table 6 shows the compensation coefficients of the system considering the geomagnetic gradient (i.e., the method proposed in this paper). A comparison of the two tables shows that most of the parameters were very close.  The level flight path is shown in Figure 10. The level flight was about 30 km, the aim of which was to verify the best performance of aeromagnetic compensation.

Data Analysis
The compensation parameters were solved according to the process shown in Figure  11. The geomagnetic gradient was considered in the process of solving the compensation coefficients. The raw data were sampled at 160 Hz, then filtered by a lowpass filter with a cut-off frequency of 20 Hz. The geomagnetic field column He was computed with Equation (16), and the information matrix A was constituted by Equation (6). Finally, the compensation parameters were solved by Equation (9). After the calibration flight, the effectiveness of the proposed method was examined based on the test flight data. The compensation coefficients could be obtained without taking into consideration the geomagnetic gradient (i.e., the traditional method), as shown in Table 5. Table 6 shows the compensation coefficients of the system considering the geomagnetic gradient (i.e., the method proposed in this paper). A comparison of the two tables shows that most of the parameters were very close.  The three-axis magnetometer data in the calibration flight are shown in Figure 12. It is obvious that there were four periods of time for the compensation maneuvers. At each stage, the pitch, yaw, and roll were executed in five cycles.  The three-axis magnetometer data in the calibration flight are shown in Figure 12. It is obvious that there were four periods of time for the compensation maneuvers. At each stage, the pitch, yaw, and roll were executed in five cycles.  Figure 13 illustrates a typical performance of compensation in four orthogonal headings. The blue lines are the uncompensated data, and the red lines are the compensated data with the proposed method. The comparison of the compensation results between the traditional model, the Ge model [28], and the proposed model are shown in Figure 14, Figure 12. The three-axis magnetometer data in the calibration flight (after a bandpass filter, offset for clarity). Figure 13 illustrates a typical performance of compensation in four orthogonal headings. The blue lines are the uncompensated data, and the red lines are the compensated data with the proposed method. The comparison of the compensation results between the traditional model, the Ge model [28], and the proposed model are shown in Figure 14, with the uncompensated signals in blue, the compensated signals of the traditional method in green, the compensated signals of the Ge method in black, and the compensated signals of the proposed method in red.
Appl. Sci. 2022, 12, 1490 13 of 16 Figure 13. The compensation results of calibration flight using coefficients of Table 6 (after a bandpass filter, offset for clarity). Five sets of data were compensated respectively. These include four sets of calibration flight data, CF1, CF2, CF3, and CF4, and one set of level flight data, marked as LF1. The compensation results obtained using the traditional method, the Ge method, and the proposed method are listed in Table 7.    Figure 13. The compensation results of calibration flight using coefficients of Table 6 (after a bandpass filter, offset for clarity). Five sets of data were compensated respectively. These include four sets of calibration flight data, CF1, CF2, CF3, and CF4, and one set of level flight data, marked as LF1. The compensation results obtained using the traditional method, the Ge method, and the proposed method are listed in Table 7.  Five sets of data were compensated respectively. These include four sets of calibration flight data, CF1, CF2, CF3, and CF4, and one set of level flight data, marked as LF1. The compensation results obtained using the traditional method, the Ge method, and the proposed method are listed in Table 7. Based on the five sets of data shown in Table 7, it was concluded that the performance of the proposed method was better than that of the traditional method and the Ge method. For the four groups of calibrated flight data, compared to the traditional method, the improvement in IR after considering the geomagnetic gradient was between 152% and 224%, which is quite significant. The improvement in IR for the level flight data reached a significant 329%. As shown in Figure 15, the peak-to-peak value of the residual noise after compensation was 19.4 pT, which proves that the proposed method can well compensate the interference field of the aircraft and improve the ability of the system to detect weak magnetic anomalies. Based on the five sets of data shown in Table 7, it was concluded that the performance of the proposed method was better than that of the traditional method and the Ge method. For the four groups of calibrated flight data, compared to the traditional method, the improvement in IR after considering the geomagnetic gradient was between 152% and 224%, which is quite significant. The improvement in IR for the level flight data reached a significant 329%. As shown in Figure 15, the peak-to-peak value of the residual noise after compensation was 19.4 pT, which proves that the proposed method can well compensate the interference field of the aircraft and improve the ability of the system to detect weak magnetic anomalies.

Conclusions
The traditional aeromagnetic compensation method does not take into consideration the influence of the geomagnetic gradient. The interference field obtained after bandpass filtering still contains the information of the geomagnetic field, which results in errors during the estimation of compensation parameters. In this paper, an improved compensation method that is robust to the geomagnetic gradient is presented. By introducing the IGRF model, the geomagnetic gradient is modeled as a function of longitude, latitude, and flight height. This significantly improved the accuracy and the stability of the estimated compensation parameters. The experiment results show that the proposed method greatly improved the compensation effect. All the results demonstrate that the proposed method could significantly improve the effect of aeromagnetic compensation and the weak magnetic anomaly detection. In the future, more aeromagnetic survey tests will be carried out to verify the effectiveness of the proposed method.

Conclusions
The traditional aeromagnetic compensation method does not take into consideration the influence of the geomagnetic gradient. The interference field obtained after bandpass filtering still contains the information of the geomagnetic field, which results in errors during the estimation of compensation parameters. In this paper, an improved compensation method that is robust to the geomagnetic gradient is presented. By introducing the IGRF model, the geomagnetic gradient is modeled as a function of longitude, latitude, and flight height. This significantly improved the accuracy and the stability of the estimated compensation parameters. The experiment results show that the proposed method greatly improved the compensation effect. All the results demonstrate that the proposed method could significantly improve the effect of aeromagnetic compensation and the weak magnetic anomaly detection. In the future, more aeromagnetic survey tests will be carried out to verify the effectiveness of the proposed method.