Compensation Method for Diurnal Variation in Three-Component Magnetic Survey

: Considering that diurnal variation interferes with three-component magnetic surveys, which inevitably a ﬀ ects survey accuracy, exploring an interference compensation method of high-precision is particularly desirable. In this paper, a compensation method for diurnal variation is proposed, the procedure of which involves calibrating the magnetometer error and the misalignment error between magnetometer and non-magnetic theodolite. Meanwhile, the theodolite is used to adjust the attitude of the magnetometer to unify the observed diurnal variation into the geographic coordinate system. Thereafter, the feasibility and validity of the proposed method were veriﬁed by ﬁeld experiments. The experimental results show that the average error of each component between the observed value of the proposed method and that of Changchun Geomagnetic station is less than 1.2 nT, which indicates that the proposed method achieves high observation accuracy. The proposed method can make up for the deﬁciency that traditional methods cannot meet the requirements of high accuracy diurnal variation compensation. With this method, it is possible for us to set up temporary diurnal variation observation station in areas with complex topography and harsh environment to assist aeromagnetic three-component survey.


Introduction
Compared with the scalar magnetic survey, three-component magnetic survey can obtain richer magnetic information, facilitate making quantitative interpretation of magnetic anomaly and enhance the resolution of magnetic target positioning, which plays an important role in geological survey, mineral exploration and earth science research [1][2][3][4][5]. However, the magnetic survey accuracy is inevitably decreased by diurnal variation, thus making it of great importance to develop an effective method for interference compensation.
The magnetic field measured by tri-axial magnetometer can be modeled as the sum of three components Where h m denotes the measured data, h g denotes the geomagnetic field, h i denotes the magnetic interference field, and h o denotes the diurnal variation interference.
In practice of three-component magnetic survey, the first component is considered as a valuable element while the other two components are as magnetic interference. The magnetic interference field 1.
Scaling error. Scaling error denotes the difference in sensitivity of each axis due to the different characteristics of the internal electronic devices. The scaling error matrix k s f can be modeled as where k 1 , k 2 , k 3 denote the axial scaling errors of the magnetometer respectively. 2.
Offset error. Offset error denotes the deviation between magnetometer's output and true value, which can be modeled as where h x b , h y b , h z b denote the axial offset errors of the magnetometer respectively. 3.
Non-orthogonality error. The orthogonality between three axes of the magnetometer cannot be guaranteed due to manufacturing precision limitations, thus resulting in non-orthogonal error, as illustrated in Figure 1. As is shown, o − xyz denotes the ideal tri-axial magnetometer while o − x y z denotes the non-orthogonal magnetometer.
The non-orthogonality error matrix k nor can be modeled as Appl. Sci. 2020, 10 = + h k k h h (4) where m h denotes the measured data of the magnetometer, g h denotes the true geomagnetic field.
According to Equation (4), the calibration model of magnetometer error can be expressed as where We can get In areas with stable magnetic field, the magnitude of the geomagnetic field remains constant within a short period, so as the magnetometer rotates, its measuring locus will form a sphere with a radius equal to the magnitude of the local geomagnetic field [11]. However, due to the magnetometer error, the sphere is distorted into an ellipsoid [12].
The equation of quadric surface is  The comprehensive mathematical model for output error of the magnetometer interfered by different error sources can be expressed as where h m denotes the measured data of the magnetometer, h g denotes the true geomagnetic field. According to Equation (4), the calibration model of magnetometer error can be expressed as In areas with stable magnetic field, the magnitude of the geomagnetic field remains constant within a short period, so as the magnetometer rotates, its measuring locus will form a sphere with a radius equal to the magnitude of the local geomagnetic field [11]. However, due to the magnetometer error, the sphere is distorted into an ellipsoid [12].
The equation of quadric surface is to the quadric surface F(ρ, σ) = 0. Therefore, the parameter ρ can be solved by fitting based on the minimum value of F(ρ, v). That is, Appl. Sci. 2020, 10, 986 4 of 13 The quadric surface has multiple forms such as ellipsoid, cylinder, parabola, etc. In the mathematical sense, to ensure that the fitted quadric surface is an ellipsoid, it is required to add the following constraints in the fitting process.
After obtaining the parameter ρ, the Equation F(ρ, v) = 0 can be written as the following matrix form Then, convert Equation (10) where Comparing Equation (7) and Equation (11), we can get Here, as A e is a positive definite matrix, Cholesky decomposition can be applied Then, we can get

Misalignment Error Calibration
In this method, a high accuracy non-magnetic theodolite is used to precisely adjust the attitude of the magnetometer, so as to unify the observed diurnal variation into the geographic coordinate system. For this purpose, we need to calibrate misalignment error between the magnetometer and the theodolite. Figure 2 illustrates misalignment error between the magnetometer and the theodolite. As is shown, we define two coordinate systems: the magnetometer coordinate system (o m − x m y m z m ) and the theodolite coordinate system (o t − x t y t z t ), in which x m , y m , z m and x t , y t , z t denote three axes of the magnetometer and the theodolite respectively.
According to the Euler-angle rotation method [13], three rotations can re-orient an object in any direction. This method can be applied to calibrate the misalignment error. Here, Appl. Sci. 2020, 10, 986

of 13
The process can be described by the following rotation matrices.
Appl. Sci. 2020, 10 x 5 of 13 where , , α β γ denote the misalignment error angles. With the known error angles, the transformation of measured data from the magnetometer coordinate system to the theodolite coordinate system is given by In this paper, we adopted an approach to identify the three misalignment error angles by rotation. When rotating the magnetometer around one axis of the theodolite, the projection of geomagnetic field on the axis of the magnetometer aligned with the rotation axis is constant after calibrating misalignment error, which serves as the basis for calibration method.
According to Equation (20), the measured data are substituted into the following equation to calibrate the misalignment error.
where m h denotes the measured data of the magnetometer rotating around one axis of the theodolite. The final rotation matrix C t m can be found by multiplying these together.
With the known error angles, the transformation of measured data from the magnetometer coordinate system to the theodolite coordinate system is given by In this paper, we adopted an approach to identify the three misalignment error angles by rotation. When rotating the magnetometer around one axis of the theodolite, the projection of geomagnetic field on the axis of the magnetometer aligned with the rotation axis is constant after calibrating misalignment error, which serves as the basis for calibration method.
According to Equation (20), the measured data are substituted into the following equation to calibrate the misalignment error.
where h m denotes the measured data of the magnetometer rotating around one axis of the theodolite.
Appl. Sci. 2020, 10, 986 As analyzed above, v x is constant during the rotation around x t -axis. Similarly, v y and v z remain unchanged as well when rotating around y t -axis and z t -axis respectively. Thus, we can get The objective function is defined as This optimization design aims to estimate three misalignment error angles by minimizing the objective function. The optimization process solves this nonlinear system with multiple objectives using the least square method [14].
Rotating the magnetometer around one axis of the theodolite can only align this axis to the magnetometer's axis of the same direction, which means rotating the magnetometer around y t -axis of the theodolite only ensures the alignment of y t -axis to y m -axis. Therefore, it is required to rotate the magnetometer around at least two axes of the theodolite respectively so as to calibrate the misalignment error. Meanwhile, it should be noted that rotation in pitch (rotating around y t -axis) of the magnetometer ought to be avoided after calibrating the misalignment error, otherwise it will invalidate the obtained misalignment error angles (α, β, γ).

Alignment to North
The geographic coordinate system with three axes pointing to north, east and up respectively is taken as the datum coordinate system for magnetic survey. Here, the north direction is the geographic North Pole, also known as true north. After finishing the above two steps of calibration, we will discuss how to align the axis of magnetometer to true north in this part.
As illustrated in Figure 3, θ is the angle between the line of two sites (A and B) of distance d and true north direction. The challenge of aligning to true north direction is to obtain the angle θ accurately. Here, the angle θ can be calculated by the following steps. Firstly, the latitude and longitude of the selected sites A and site B are measured using differential GPS. Then, the latitude L and longitude B are converted into x and y coordinate values of plane coordinates by means of coordinate transformation to obtain the relative position of sites A and B [15]. The transformation of the coordinate system is as follows.
with ' a ' denoting the semi-major axis of the earth, ' F ' denoting the flattening factor of the earth. The prime vertical radius of curvature N is where the first eccentricity of earth e is   We can substitute the coordinate values of sites A and B into Equation (30) to calculate the angle θ . Theoretically, the farther the distance between sites A and B, the higher the accuracy of angle θ is. A distance between the two sites no less than 100 m is acceptable.
where the first eccentricity of earth e is The transformation equation is given by with 'L 0 ' denoting the longitude of central meridian. Where t = tan B, η = e 2 1−e 2 , y o =500,000. We can substitute the coordinate values of sites A and B into Equation (30) to calculate the angle θ. Theoretically, the farther the distance between sites A and B, the higher the accuracy of angle θ is. A distance between the two sites no less than 100 m is acceptable.
where (x A , y A ) and (x B , y B ) are coordinate values of the sites A and B respectively. Before aligning the axis of magnetometer to true north direction using the obtained angle θ, the following preparations are quite necessary: (1) place the theodolite on the site B using the optical plummet; (2) level the theodolite by adjusting the leveling foot screw; (3) align the reticle of objective lens to site A by adjusting the theodolite. After making the above preparations, the coming procedure is to rotate the theodolite to align the axis of magnetometer to true north. It should be noted that the relative position of sites A and B will affect the rotation angle of the theodolite. There are four different situations for the relative position of sites A and B, which can be expressed in quadrants as: NE-quadrant, NW-quadrant, SW-quadrant, and SE-quadrant, as illustrated in Figure 4. The rotation angle = θ i (i = 1, 2) when the relative position is in NE-quadrant and NW-quadrant while rotation angle = 180 0 − θ i (i = 3, 4) in SW-quadrant and SE-quadrant.
is to rotate the theodolite to align the axis of magnetometer to true north. It should be noted that the relative position of sites A and B will affect the rotation angle of the theodolite. There are four different situations for the relative position of sites A and B, which can be expressed in quadrants as: NE-quadrant, NW-quadrant, SW-quadrant, and SE-quadrant, as illustrated in Figure 4

Experimental Setup
In order to verify the feasibility and validity of the proposed method, a diurnal variation observation system consisting of a tri-axial magnetometer (MAG-03MS100, Bartington, Witney, Britain), a non-magnetic theodolite (TDJ2E-NM, Boif, Beijing, China) and a differential GPS (IMU-IGM-S1, Novatel, Calgary, AB, Canada) was built in the experiment. The performance specifications of the experimental setup are shown in Table 1.

Experimental Setup
In order to verify the feasibility and validity of the proposed method, a diurnal variation observation system consisting of a tri-axial magnetometer (MAG-03MS100, Bartington, Witney, Britain), a non-magnetic theodolite (TDJ2E-NM, Boif, Beijing, China) and a differential GPS (IMU-IGM-S1, Novatel, Calgary, AB, Canada) was built in the experiment. The performance specifications of the experimental setup are shown in Table 1. As described previously, the magnetometer error, misalignment error, and alignment to the north are calibrated successively. The relevant experiment was carried out in an area within 100 m far from Changchun Geomagnetic station. There is no magnetic interference around the selected area. In addition, the six-channel spectrum analyzer Spectramag-6 supporting the use of Mag-03MS100 magnetometer is adopted as the data acquisition unit.

Calibration Results
The total field intensity of experimental site measured by GSM-19 Overhauser magnetometer (GEM Systems, Markham, ON, Canada) is 55,046.65 nT, which is used as a reference for calibrating the magnetometer error. During the process of data collection, the magnetometer is rotated to capture samples with different attitudes. From Figure 5a, it can be seen that interfered by the magnetometer error, the measurement locus appears as an ellipsoid. With the fitted ellipsoid, we are capable of obtaining the magnetometer error calibration parameters k and h b . Figure 5b illustrates the calibration result, from which we can see that the maximum residual of total field is reduced from 105 nT to less than 3 nT, indicating a significant decrease in the measurement error of the magnetometer.
As described previously, the magnetometer error, misalignment error, and alignment to the north are calibrated successively. The relevant experiment was carried out in an area within 100 m far from Changchun Geomagnetic station. There is no magnetic interference around the selected area. In addition, the six-channel spectrum analyzer Spectramag-6 supporting the use of Mag-03MS100 magnetometer is adopted as the data acquisition unit.

Calibration Results
The total field intensity of experimental site measured by GSM-19 Overhauser magnetometer (GEM Systems, Markham, ON, Canada) is 55,046.65 nT, which is used as a reference for calibrating the magnetometer error. During the process of data collection, the magnetometer is rotated to capture samples with different attitudes. From Figure 5a, it can be seen that interfered by the magnetometer error, the measurement locus appears as an ellipsoid. With the fitted ellipsoid, we are capable of obtaining the magnetometer error calibration parameters k and b h . Figure 5b illustrates the calibration result, from which we can see that the maximum residual of total field is reduced from 105 nT to less than 3 nT, indicating a significant decrease in the measurement error of the magnetometer.
(a) (b) In the misalignment error calibration experiment, rotations of the magnetometer around t yaxis and t z -axis are performed separately to obtain experimental data. As mentioned above, the pitch angle of the magnetometer ought to avoid being changed after calibrating misalignment error. Hence, the rotation sequence is first around t y -axis and then t z -axis. Figure 6 illustrates the misalignment error calibration result. It can be seen that, interfered by the misalignment error during rotation, the z /nT In the misalignment error calibration experiment, rotations of the magnetometer around y t -axis and z t -axis are performed separately to obtain experimental data. As mentioned above, the pitch angle of the magnetometer ought to avoid being changed after calibrating misalignment error. Hence, the rotation sequence is first around y t -axis and then z t -axis. Figure 6 illustrates the misalignment error calibration result. It can be seen that, interfered by the misalignment error during rotation, the fluctuation of y-component exceeds 981.1 nT while that of z-component exceeds 269.2 nT. Meanwhile, we can see that the fluctuations have been reduced to less than 7.4 nT and 5.9 nT respectively after calibration, which indicates the misalignment error is decreased to an acceptable range. Meanwhile, we can see that the fluctuations have been reduced to less than 7.4 nT and 5.9 nT respectively after calibration, which indicates the misalignment error is decreased to an acceptable range. Prior to aligning the axis of the magnetometer to the true north direction, the theodolite is required to be leveled. With the help of a high-precision plate level, the horizontal error can be restricted within 20″ after leveling. In the procedure of aligning to the north, the latitude and longitude of the selected sites A and B are measured as (44.08735182° N, 124.90305218° E) and (44.0875215302° N, 124.90305218° E), as illustrated in Figure 7. The angle θ is calculated to be 83.49442° using the above positional parameters. Meanwhile, the relative position of the two sites is in SW-quadrant, which means the theodolite requires to be rotated 96.50558° in the north direction to align the t x -axis of magnetometer to the true north.
After the above steps of calibration, it can be assumed that the magnetometer error and misalignment error have been eliminated and the t x -axis of the magnetometer is aligned to the true north direction.  Prior to aligning the axis of the magnetometer to the true north direction, the theodolite is required to be leveled. With the help of a high-precision plate level, the horizontal error can be restricted within 20" after leveling. In the procedure of aligning to the north, the latitude and longitude of the selected sites A and B are measured as (44.08735182 • N, 124.90305218 • E) and (44.0875215302 • N, 124.90305218 • E), as illustrated in Figure 7. The angle θ is calculated to be 83.49442 • using the above positional parameters. Meanwhile, the relative position of the two sites is in SW-quadrant, which means the theodolite requires to be rotated 96.50558 • in the north direction to align the x t -axis of magnetometer to the true north. respectively after calibration, which indicates the misalignment error is decreased to an acceptable range. Prior to aligning the axis of the magnetometer to the true north direction, the theodolite is required to be leveled. With the help of a high-precision plate level, the horizontal error can be restricted within 20″ after leveling. In the procedure of aligning to the north, the latitude and longitude of the selected sites A and B are measured as (44.08735182° N, 124.90305218° E) and (44.0875215302° N, 124.90305218° E), as illustrated in Figure 7. The angle θ is calculated to be 83.49442° using the above positional parameters. Meanwhile, the relative position of the two sites is in SW-quadrant, which means the theodolite requires to be rotated 96.50558° in the north direction to align the t x -axis of magnetometer to the true north.
After the above steps of calibration, it can be assumed that the magnetometer error and misalignment error have been eliminated and the t x -axis of the magnetometer is aligned to the true north direction.  After the above steps of calibration, it can be assumed that the magnetometer error and misalignment error have been eliminated and the x t -axis of the magnetometer is aligned to the true north direction.

Diurnal Variation Observation
In order to evaluate the performance of the proposed method, a field experiment is conducted subsequently. The diurnal variation of 24 h was observed with a sampling frequency of 1 Hz in the experiment. Meanwhile, the diurnal variation obtained by the Changchun Geomagnetic station was used as the reference for evaluating the observation accuracy of the proposed method. Figure 8 illustrates the diurnal variation observations. It can be intuitively seen that the trend of observations of the proposed method is highly consistent with the reference. Besides, two indexes of Pearson correlation coefficient p [16] and average error e (Equation (28)) are introduced in the paper to evaluate the performance of the proposed method more objectively, as it is shown in Table 2. The larger p is, the higher the correlation is and the smaller e is, the higher the accuracy of observation is. Commonly, 0.8 ≤ p < 1 denotes high correlation. The fact that p of each component is greater than 0.8 and e is less than 1.2 nT indicates high observation accuracy of the proposed method.
where o denotes observed value and o denotes the average value of o, r denotes reference value and r denotes the average value of r.

Diurnal Variation Observation
In order to evaluate the performance of the proposed method, a field experiment is conducted subsequently. The diurnal variation of 24 h was observed with a sampling frequency of 1 Hz in the experiment. Meanwhile, the diurnal variation obtained by the Changchun Geomagnetic station was used as the reference for evaluating the observation accuracy of the proposed method. Figure 8 illustrates the diurnal variation observations. It can be intuitively seen that the trend of observations of the proposed method is highly consistent with the reference. Besides, two indexes of Pearson correlation coefficient p [16] and average error e (Equation (28)) are introduced in the paper to evaluate the performance of the proposed method more objectively, as it is shown in Table 2. The larger p is, the higher the correlation is and the smaller e is, the higher the accuracy of observation is. Commonly, 0.

Conclusions
Aiming at the reduction of the diurnal variation in three-component magnetic survey, a compensation method is proposed in the paper. With this method, the diurnal variation in the geographic coordinate system is accurately observed, based on which, the interference can be effectively compensated using the obtained magnetic variation data. The feasibility and validity of the method are verified by field experiment, the results of which indicate that the proposed method can meet the requirement of three-component magnetic survey for diurnal variation compensation accuracy. Compared with the traditional methods, the proposed method overcomes the deficiency of immovability of fixed geomagnetic stations and low compensation accuracy of modeling approach. With the advantages of simple operation and portability of equipment used, the proposed method can be used to assist aeromagnetic three-component survey in areas with special environment such as deserts and forests by setting up the temporary observation station.