Accounting for Geometric Anisotropy in Sparse Magnetic Data Using a Modified Interpolation Algorithm

Abstract

The construction of a high-precision geomagnetic map is a prerequisite for geomagnetic navigation and magnetic target-detection technology. The Kriging interpolation algorithm makes use of the variogram to perform linear unbiased and optimal estimation of unknown sample points. It has strong spatial autocorrelation and is one of the important methods for geomagnetic map construction. However, in a region with a complex geomagnetic field, the sparse geomagnetic survey lines make the ratio of line-spacing resolution to in-line resolution larger, and the survey line direction differs from the geomagnetic trend, which leads to a serious effect of geometric anisotropy and thus, reduces the interpolation accuracy of the geomagnetic maps. Therefore, this paper focuses on the problem of geometric anisotropy in the process of constructing a geomagnetic map with sparse data, analyzes the influence of sparse data on geometric anisotropy, deduces the formula of geometric anisotropy correction, and proposes a modified interpolation algorithm accounting for geometric anisotropy correction of variogram for sparse geomagnetic data. The results of several sets of simulations and measured data show that the proposed method has higher interpolation accuracy than the conventional spherical variogram model in a region where the geomagnetic anomaly gradient changes sharply, which provides an effective way to build a high-precision magnetic map of the complex geomagnetic field under the condition of sparse survey lines.

1. Introduction

As a basic work of geomagnetic matching navigation and magnetic anomaly detection, the geomagnetic map is of great significance to improve navigation accuracy and target-detection ability. Geomagnetic mapping can be obtained by means of aeromagnetic survey, etc. According to relevant measurement specifications [1,2,3,4,5], the magnetic data obtained by aeromagnetic survey is characterized by dense data in the direction of survey lines and sparse data in the direction of line-spacing. Such data cannot be directly used for geomagnetic navigation and magnetic target detection, and it is necessary to interpolate the data to establish high-resolution grid data. The Kriging interpolation algorithm is one of the best methods. Kriging uses a variogram to perform linear unbiased and optimal estimation of unknown sample points, and it has strong spatial autocorrelation, which can better represent the structural component of the geomagnetic field and accurately reflect the geomagnetic information in a local region [6]. Hao et al. [7] studied the basic principle of ordinary Kriging interpolation and found that the accuracy of geomagnetic grid data obtained using it is higher than that obtained by inverse distance weighting interpolation, natural neighbor interpolation and other methods. Yang et al. [8] used cross-validation to verify that the universal Kriging method has better interpolation effect than the ordinary Kriging method. Cai et al. [9] applied the co-Kriging interpolation algorithm to indoor navigation and positioning based on the geomagnetic field. The sample was divided into two parts for ordinary Kriging interpolation processing and then fusion processing, which has the characteristics of small root-mean-square error and good interpolation effect. At the same time, with the development of deep learning technology, Zhao et al. [10] combined neural networks with a Kriging interpolation algorithm and proposed a geomagnetic reference map construction method based on the combination of Kriging and a back propagation (BP) neural network. The Kriging interpolation error in a large gradient region is greater than that in a BP neural network, and the interpolation error in a small gradient region is smaller than that in a BP neural network. In the large gradient region, the predicted value of the BP neural network is used instead of the Kriging interpolation result to improve the interpolation accuracy.

However, the sparsity of aeromagnetic data and the direction of survey line intensify the influence of geometric anisotropy on the accuracy of geomagnetic map construction, which shows that the geomagnetic detail with high resolution in-line cannot be distinguished with low resolution in the line-spacing direction. Wackernagel et al. [11], Sun [12] and Mcbratney et al. [13] point out that anisotropy is absolute and isotropy is relative. However, geometric anisotropy is often ignored in the actual processing, so the accuracy of interpolation is reduced by processing according to isotropy. Therefore, this paper analyzes the geometric anisotropy problem of sparse data, deduces the geometric anisotropy processing formula, and proposes a modified interpolation algorithm accounting for geometric anisotropy correction of variogram for sparse geomagnetic data. Based on the simulated data and the measured data, the proposed algorithm is used to construct high-precision geomagnetic maps respectively, and is compared with the geomagnetic maps constructed using a spherical variogram model. The results show that the interpolation accuracy of this method has obvious advantages, which greatly reduces the maximum absolute error, average absolute error and root-mean-square error, and is suitable for constructing local geomagnetic reference maps with high precision and high resolution.

2. Influence Analysis of Geometric Anisotropy of Variogram

2.1. Principles of Kriging Interpolation Algorithm

The Kriging interpolation algorithm is a spatial interpolation method based on geostatistics, which makes a minimum variance estimation and linear unbiased method for unknown data of similar characteristics in a region, according to certain characteristic data of several information samples.

2.1.1. Regionalized Variable

The magnetic data provide information about regionalized variables, which are simply functions $z(x,y)$ (of a spatial or time continuum) whose behavior we would like to characterize. In the probabilistic model, a regionalized variable $z(x,y)$ is considered to be a realization of a random function $Z(x,y)$ (the infinite family of random variables constructed at all points x. of a spatial domain). The advantage of this approach is that we shall only try to characterize the random function $Z(x,y)$ and not a particular regionalized variable. In such a setting, the data values are samples from one particular realization $z(x,y)$ of $Z(x,y)$ [11].

2.1.2. Principle of Interpolation

The basic idea of the Kriging algorithm is that a certain weight is assigned to each sample value, according to the spatial correlation between sample points, especially the spatial position, and the unknowns in the estimation region that were estimated using the weighted average method [14].

As shown in Figure 1, $Z_i(i=1,2,\dots ,n)$(the red points) is a set of measured magnetic field sample data points with coordinate position $(x_i,y_i)$, and the estimated value ${\widehat{Z}}_0$(the green point) of a prediction point with coordinate position $(x_0,y_0)$ can be expressed by linear combination of sampled magnetic field samples:

$${\widehat{Z}}_0={\displaystyle \sum }_{i=1}^n{\lambda }_i{Z}_i$$

(1)

where ${\lambda }_i$ is the weight of the measured point $Z_i$, ${\widehat{Z}}_0$ is called the Kriging estimation of $Z_0$. As long as the weight coefficient ${\lambda }_i$ is determined, the forecast point estimation ${\widehat{Z}}_0$ can be obtained. Kriging interpolation estimation is used to calculate the weight ${\lambda }_i$ under the condition that the estimate ${\widehat{Z}}_0$ is unbiased and the estimation variance is minimum.

2.1.3. Variogram

The variogram is the core of the Kriging interpolation algorithm and is defined as half of the variance of the difference of the regionalized variable $Z(x)$ at two points $x,x+l$, and the formula is:

$$\begin{array}{ll}\gamma (l)=& \frac{1}{2}Var\left[Z(x)-Z(x+l)\right]\\ & =\frac{1}{2}\left\{E{[Z(x)-Z(x+l)]}^2\right.\left.-{(E[Z(x)-Z(x+l)])}^2\right\}\\ & =\frac{1}{2}E{[Z(x)-Z(x+l)]}^2\end{array}$$

(2)

where $\gamma (l)$ is the variogram, and $l$ is the lag distance. There are many variogram models, such as the spherical, which is more commonly used; its general form is as follows:

$$\gamma (l)=\left\{\begin{array}{cc}{C}_0,& l=0\\ C_0+C(\frac{3l}{2a}-\frac{l^3}{2a^3}),& 0<l\le a\\ C_0+C& l>a\end{array}\right.$$

(3)

where $C_0$ is the nugget, which is the value of the variogram, when the lag $l$ is 0; $C_0+C$ is the sill, which is defined as the value reached when the variogram is in a stable state with the increase of lag $l$; $C$ is the arch height, which is defined as the difference between sill and nugget; $a$ is the range, which is defined as the size of the autocorrelation range of regionalized variables.

2.1.4. Experimental Variogram

The experimental variogram ${\gamma }^{\ast }(l)$ is an estimate of the variogram $\gamma (l)$ constructed from the measured data. We can divide every pair of data: $\left\{Z(x_i),Z(x_i+l)\right\}$ $(i=1,2,\dots ,N(l))$ is respectively regarded as a one-time implementation of $\left\{Z(x_i),Z(x_i+l)\right\}$, $N(l)$ represents the groups of the data pair divided by the vector $l$, so the formula of the experimental variogram is obtained:

$${\gamma }^{\ast }(l)=\frac{1}{2N(l)}{\displaystyle \sum }_{i=1}^{N(l)}{[Z(x_i)-Z(x_i+l)]}^2$$

(4)

For different lags $l$, the values of ${\gamma }^{\ast }(l)$ can be calculated according to Equation (4), and points $[l_i,\gamma (l_i)]$ can be marked in the graph of experimental variogram, which can directly reflect the relevant situation of measurement data [11].

2.2. Analysis of Geometric Anisotropy of Experimental Variogram

The constructed geomagnetic map is a two-dimensional plane with the height $h$ as the constant value, $u$ and $v$ as the coordinate axes, and its surveying simulation is shown in Figure 2. The $u$ and $v$ axes, which are perpendicular to each other, represent the direction in the flight-line and the direction of line-spacing during aeromagnetic measurement, and the circles with black dot are the sampling points. The magnetic data value is a regionalized variable, and the selected experimental variogram (such as spherical variogram) diagram can be drawn by sample points in Figure 3. $a_u$ and $a_v$ are the variables in different directions of the $u$ and $v$ axis, and $K=a_u/a_v$ is the anisotropy ratio. When $K=1$, the direction of $u$ and $v$ axes is isotropic; when $K\ne 1$, $a_u<a_v$ (the blue solid line in Figure 2) or $a_u>a_v$ indicates the direction of $u$ and $v$ axes is anisotropic. Its physical meaning is: the average degree of variation between two points with a distance of $l$ in the direction of $u$ is equal to that with a distance of $Kl$ in the direction of $v$, that is, ${\gamma }_u(l)={\gamma }_v(Kl)$.

Using this variogram model, the diagram of direction–range is shown in Figure 4a. When $K=1$, the direction of $u$ and $v$ axes is isotropic, the graph of direction–range is approximately a circle with range $a_u$ (or $a_v$) as the radius, which is shown by the dotted line in Figure 4a. When $K\ne 1$, the direction of $u$ and $v$ axes is anisotropic, the graph of direction–range is approximately an ellipse with range $a_u$ as the long radius and $a_v$ as the short radius, which is shown by the solid line in Figure 4a. This can be understood as meaning that the characteristic direction of the geomagnetic data is inconsistent with the mapping direction. If the azimuth angle between the geomagnetic variation trend and the measured $u,v$ axis is $\theta$, then the variogram diagram and the direction–range diagram will change due to the introduction of $\theta$. The converted variogram diagram is shown by the red line in Figure 3, and the converted direction–range diagram is shown by the solid line in Figure 4b. The red and blue lines in Figure 3 correspond to the red and blue lines in Figure 4b, indicating the transformation of the coordinate system from Coordinate $(u,v)$ to Coordinate $(u^{\prime },v^{\prime })$.

Therefore, the results of geomagnetic data interpolation are not only related to the length of lag, but also to the direction of lag. If the trend of geomagnetic change is the same as the change of lag distance in all directions, then it is isotropic. If the geomagnetic trend is not the same as the change of the lag distance in each direction, it shows geometric anisotropy. This is caused by the different evolution of underlying physical processes in space [15].

For sparse geomagnetic data, using an isotropic interpolation method to process anisotropic geomagnetic data has a great impact on the interpolation accuracy. For further explanation, this paper draws the interpolation processing diagram using the known sampling points, as shown in Figure 5. The figure shows the variable range and orientation relationship between the center point and other points in the effective region (the big pink circle region); the black dots represent the sampling points; the sampling mode along the $u$ and $v$ directions is to simulate the sampling mode between the line-spacing direction ($u$ direction) and the flight-line direction ($v$ direction) during aeromagnetic survey in Figure 2; and the sampling points along the  direction are twice the number of sampling points along the $u$ direction to simulate the sparsity of geomagnetic data. Although two times is a small percentage, it is just one way of analyzing the problem.

It can be seen from Figure 5 that: (a) through the length of the arrows of different colors, the size of the interval between the center point and the sampling points in each direction changes with the change of the direction angle; and (b) by setting the effective region between two points (the large, medium and small regions of the pink circle), the number of sampling points in each direction around the central point decreases with the reduction of the effective region, and some directions even have no sample points.

It can be determined from Figure 5 that:

(1)

If the geomagnetic trend is distributed along the direction $v$ (orange arrows), the number of sampling data is sufficient and the interval is small, thus the resolution is high and the magnetic trend in this direction can be better represented;

(2)

If the geomagnetic trend is distributed along the direction $u$ (red arrows), the number of sampling data and interval are significantly lower than that along the direction $v$, which can properly reflect the geomagnetic trend in this direction. If the sampling interval in this direction is increased, and the data become more and more sparse, its resolution will be reduced and the interpolation accuracy will be reduced, too;

(3)

If the geomagnetic trend is distributed along the blue line, the number of sampling points in this direction is extremely reduced (the number of blue arrows), and the sampling interval is significantly larger than the above two methods, which makes the details of geomagnetic changes in this direction easy to be ignored, thus reducing the accuracy of interpolation. At this time, in order to increase the number of samples, the area in the azimuth direction is often expanded by setting the tolerance $\alpha$ in the azimuth direction $\theta$ and the lag distance tolerance, as shown in Figure 6 for the blue direction expansion in Figure 5, where the azimuth angle $\theta =56°$ and the azimuth tolerance $\alpha =15°$.

The Figure 5 and Figure 6 show the method of sample point selection before the calculation of variogram parameters in the interpolation processing of geomagnetic data, which directly affects the calculation of the experimental variogram. As can be seen from Figure 6, during the measurement process, the recommended direction of the survey line is consistent with the trend of local geomagnetic change, and the actual sampling points are dense and the lag distance is small, which can better reflect the geomagnetic change in a small area. However, in the actual measurement process, it is difficult to match the geomagnetic trend, and the size of the interval between lines in direction $u$ is 8–10 times that in direction $v$, thus increasing the influence of geometric anisotropy on the geomagnetic construction accuracy. Therefore, this paper proposes a geomagnetic map construction algorithm based on geometric anisotropy correction of variogram, by calculating the variogram in each direction of the layout area and searching for the variogram in the direction that best fits the trend of the geomagnetic map to be processed to improve the interpolation accuracy of the geomagnetic data.

3. Construction Method of Geomagnetic Map Based on Geometric Anisotropy Correction of Variogram

The geometric anisotropic variogram is a universal problem, which has a greater impact on sparse geomagnetic data. In this section, a geomagnetic interpolation algorithm based on geometric anisotropy correction of variogram is proposed by deducing the geometric anisotropy correction formula, which provides effective support for the construction of high-precision geomagnetic maps.

3.1. Geometric Anisotropy Correction Formula

For the geometric anisotropy mentioned in Section 2.2, the relevant formulas can be deduced into the following two cases:

• The characteristic direction of the geomagnetic map is consistent with the surveying direction

At this time, the diagram of direction–range is shown by the solid line in Figure 4a. Using the sample data, the range $a_u,a_v$ in different directions $u,v$ are calculated respectively, and then updated by the following formula:

$$l=\sqrt{{l_u}^2+{(K{l}_v)}^2}$$

(5)

where anisotropy ratio $K=a_u/a_v$.

2.

The characteristic direction of the geomagnetic map is inconsistent with the surveying direction

At this time, the diagram of direction–range is shown in Figure 4b. The surveying direction (selected axis) is at the azimuth angle $\theta$ to the characteristic direction of the geomagnetic data (elliptic principal axis). The processing goes through three processes of coordinate rotation, anisotropy processing and coordinate restoration in turn. For the convenience of calculation, it is converted into matrix form and solved as follows:

$$l^{\prime }=\left[\begin{array}{c}{l^{\prime }}_u\\ {l^{\prime }}_v\end{array}\right]=Al=A\left[\begin{array}{c}{l}_u\\ l_v\end{array}\right]$$

(6)

where $A$ is the rotation matrix, $A=R_{-\theta }Q{R}_{\theta }$; $R_{\theta }$ is the coordinate system transformation matrix, $R_{\theta }=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]$; $Q$ is the anisotropic processing matrix, $Q=\left[\begin{array}{cc}1& 0\\ 0& K\end{array}\right]$; $K=a_u/a_v$; $R_{-\theta }$ is the coordinate system to restore the matrix, $R_{-\theta }=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]$ then

$$A=\left[\begin{array}{cc}{\cos }^2\theta +K{\sin }^2\theta & (1-K)\sin \theta \cos \theta \\ (1-K)\sin \theta \cos \theta & {\sin }^2\theta +K{\cos }^2\theta \end{array}\right]$$

(7)

Suppose $a_1=({\cos }^2\theta +K{\sin }^2\theta )$, $a_2=a_3=(1-K)\sin \theta \cos \theta$, $a_4={\sin }^2\theta +K{\cos }^2\theta$, then bring them into formula (7).

$$\begin{array}{ll}{l^{\prime }}^2& ={l^{\prime }}_u^2+{l^{\prime }}_v^2={(a_1{l}_u+a_2{l}_v)}^2+{(a_3{l}_u+a_4{l}_v)}^2\\ & =({a_1}^2+{a_3}^2){l_u}^2+2(a_1{a}_2+a_3{a}_4)l_u{l}_v+({a_2}^2+{a_3}^2){l_v}^2\\ & ={l_u}^2+{l_v}^2+(K^2-1){(\sin \theta l_u+\cos \theta l_v)}^2\end{array}$$

(8)

From formula (8), when $K=1$, it is geometric isotropy; when $K\ne 1$, $\theta =0°$, it is equal to formula (5), in which the characteristic direction of the geomagnetic map is consistent with the mapping direction.

3.2. Geomagnetic Interpolation Algorithm Based on Geometric Anisotropy Correction of Variogram

In the previous section, the geometric anisotropy correction formula was deduced. Based on this formula, a modified interpolation algorithm for sparse geomagnetic data is proposed in this section, and its flow chart is shown in Figure 7.

As shown in Figure 7, the algorithm mainly consists of two steps: data increase processing in the direction perpendicular to the line using the flight-line interval and construction of the magnetic map using the increased and raw magnetic data.

(1)

Data Increase Processing

According to the above sparse data survey methods, the aeromagnetic data are dense in the survey-line and sparse in the line-spacing direction, which is perpendicular to the survey-line. This step is used to increase sparse data in the line-spacing direction with the same resolution as in the survey-line, at which point a grid is constructed to reduce the effects of geometric anisotropy in both directions.

(1) According to the distribution of surveying data, the interval of the survey-line is $\alpha$ meters, the interval of the line-spacing is $\beta$ meters, and the number of the increased rows between line-spacing is calculated using equation $k=[\alpha /\beta ]-1$, in which $k$ is the number of the increased rows.

(2) Calculate the angle $\theta$ between the trend analysis of the regional geomagnetic field and the mapping direction and construct the variogram. In this paper, the brute force search method is used for processing, which mainly includes three steps: First, calculate the corresponding parameters $C_0,C,a$ of the variogram using Equation (3), and the range obtained is. Secondly, calculate the corresponding parameters $C_0,C,a$ of the variogram in each direction using Formula (3), too. The range obtained at this time is, so as to calculate the anisotropy ratios $K=a_1/a_2$ in each direction. Finally, the variogram curves of each direction and sample points are drawn in the same variogram diagram to find out which curve better fits the distribution requirements of the sample points; its corresponding azimuth angle $\theta$ is the appropriate one, and its corresponding parameters $C_0,C,a$ are the desired curve. In this paper, the initial value of the azimuth angle $\theta$ is set to 0°, the maximum value is 180°, and the step angle $\Delta \theta$ is 5° (or 6°).

(3) Based on the results of $\theta ,C_0,C,a,\mathrm{K}$ calculated in the previous step, calculate the magnetic values of the column data obtained in the first step using Kriging interpolation processing through Equations (3) and (10), so as to obtain the increased geomagnetic map data.

(2)

Construction of Geomagnetic Map

Using the increased data and known data as the sample, the data at this time fits a regular grid. Next, the Kriging method is used to fit it again, according to the required resolution, so as to construct high-precision geomagnetic map data. Specific methods can be referred to in Oliver [16].

In short, in order to construct the high-precision geomagnetic map, this paper first uses geometric anisotropy correction to increase the number of sparse data, and then uses the regular grid sample data composed of the known data and the increased data to interpolate again, according to the required resolution through the Kriging algorithm.

4. Simulation

According to Chen et al. [17], the simulation of the magnetic field is obtained using the mathematical expression of the field source with regular geometry distributed in space, such as the sphere model. This paper uses the sphere magnetic source model (referred to as the sphere combination model) to calculate, and the magnetic anomaly expression generated by a single sphere magnetic source is as follows:

$$\begin{array}{ll}\Delta T\left(x,y,z\right)& =\frac{{\mu }_0}{4\pi }\frac{m}{{\left[x_1^2+y_1^2+z_1^2\right]}^{5/2}}[(2z_1^2-x_1^2-y_1^2){\sin }^{2}I+(2x_1^2-y_1^2-z_1^2){\cos }^{2}I{\cos }^{2}D\\ & +\left(2y_1^2-x_1^2-z_1^2\right){\cos }^{2}I{\sin }^{2}D\hspace{0.33em}-3x_1{z}_1\sin 2I\cos D\hspace{0.33em}+3x_1{y}_1{\cos }^{2}I\sin 2D\\ & \left.-3y_1{z}_1\sin 2I\sin D\right]\end{array}$$

(9)

where ${\mu }_0$ represents the vacuum permeability; $D$ represents the deflection angle of the magnetization intensity; and $I$ indicates the inclination of the geomagnetic field and the inclination of magnetization intensity. $m$ represents the magnetic moment, which has a relationship with the magnetic susceptibility $J$:$m=J\times \pi R^3$, $R$ represents the radius of the sphere. $\Delta T\left(x,y,z\right)$ represents the magnetic anomaly at the position point $\left(x,y,z\right)$, the variable $\left(x_1,y_1,z_1\right)$ is calculated by the formula $x_1=x-x_0$, $y_1=y-y_0$, $z_1=z-z_0$ and $\left(x_0,y_0,z_0\right)$ represents the center coordinates of the sphere.

4.1. Simulation of Geomagnetic Field Using Two Magnetic Dipoles

In order to verify the universality of the calculation, two magnetic dipoles are used to simulate the local geomagnetic anomaly region which is shown in Figure 8a. The parameters are shown in

Table 1. Parameters of two magnetic dipoles.

No.	Position			Radius of the Sphere/Meters
	X/Meters	Y/Meters	Z/Meters
1	200	2900	−4000	400
2	7000	7500	−4400	400

. The z-axis is specified as the positive direction, the normal field of geomagnetic background is set to 48,000 nT, and the true value of the simulated geomagnetic map with a regional resolution of 50 × 50 m is generated according to the relevant parameters. The magnetic data distribution shown in Figure 8b was obtained by aeromagnetic survey, where the spacing-line interval was 1000 m and the in-line interval was 250 m.

As can be seen from Figure 9a, the trend of the geomagnetic map is about 34° upward along the X direction, by calculating the relative position of two magnetic dipoles’ points (200, 2900) and (7000, 7500). Through data increase processing, the modified algorithm selected azimuth angle 35° to be the most suitable for geometric anisotropy processing. As mentioned above, it is necessary to increase the data of four routes, as shown in Figure 9d. After processing the data with the algorithm proposed in this paper, the results are shown in Figure 9e, and the comparison results between simulation diagram and interpolation diagram are shown in Figure 9f. At the same time, to verify the effectiveness of the algorithm, we use the traditional spherical variogram to process it, and the results are shown in Figure 9b. The comparison results between the simulation graph and the interpolation graph are shown in Figure 9c.

At the same time, in order to verify the effectiveness of the algorithm, we use the maximum absolute error (ME), the mean absolute error (MAE), and the root-mean-square error (RMSE) to evaluate the accuracy. ME responses are the magnitude of the error and MAE and RMSE response distribution of the error; they are calculated as shown in Equation (10).

$$\begin{gathered} ME=\underset{i\in (1,n)}{\max }\left|z_i-{\widehat{z}}_i\right| \\ MAE=\frac{1}{n}{\displaystyle \sum _{\mathrm{i}=1}^n\left|z_i-{\widehat{z}}_i\right|} \\ RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{\mathrm{i}=1}^n{\left(z_i-{\widehat{z}}_i\right)}^2}} \end{gathered}$$

(10)

where the parameter $n$ is the number of points in Kriging interpolation, the parameter $z_i$ is the true value of the magnetic field, and the parameter ${\widehat{z}}_i$ is the magnetic value after interpolation. From the comparison (as shown in

Table 2. Comparison of Kriging interpolation errors for the simulated geomagnetic map with two magnetic dipoles using different variogram models.

Variogram Model	ME (nT)	MAE (nT)	RMSE (nT)
Spherical Variogram	40.6872	2.9604	4.7042
Improved Variogram	4.54835	0.6291	0.7985

), it can be seen that the ME, MAE, and RMSE of the improved variogram have decreased.

4.2. Simulation of Geomagnetic Field Using Multiple Magnetic Dipoles

The above simulated geomagnetic map is not representative. Therefore, six magnetic dipoles are used in this paper to simulate the local geomagnetic anomaly region. The relevant parameters are shown in

Table 3. Parameters of the six spherical combined model.

No.	Position			Radius of the Sphere/Meters
	X/Meters	Y/Meters	Z/Meters
1	1000	7000	−1000	80
2	5000	4500	−200	63
3	4500	6000	−200	63
4	4000	4500	−250	63
5	4000	5500	−200	63
6	5000	5000	−200	14

. The z-axis is specified as the positive direction, the normal field of geomagnetic background is set to 48,000 nT, and the true value of the simulated geomagnetic map with a regional resolution of 50 × 50 m is generated according to the relevant parameters, as shown in Figure 10a. In order to obtain the magnetic samples in the flight-line during the simulation of the aeromagnetic survey, the magnetic data of simulated flight along the Y direction is taken as the measured value on the constructed geomagnetic map, the in-line resolution is set to 250 m, and the flight line spacing is set to 500 m. The schematic diagram of the measurement points is shown in Figure 10b.

As mentioned above, through data increase processing, the modified algorithm-selected coordinate azimuth angle of 105° is selected for geometric anisotropy processing. The results after processing using the algorithm proposed in this paper are shown in Figure 11e, and the comparison results between the simulation graph and the interpolation graph are shown in Figure 11f. At the same time, this paper uses the traditional spherical variogram to process it, and the results are shown in Figure 11b. The comparison results between the simulation graph and the interpolation graph are shown in Figure 11d. The calculated maximum absolute error (ME), mean absolute error (MAE) and root-mean-square error (RMSE) are shown in

Table 4. Comparison of Kriging interpolation errors for the simulated geomagnetic map with six magnetic dipoles using different variogram models.

Variogram Model	ME (nT)	MAE (nT)	RMSE (nT)
Spherical Variogram	229.81627	1.6513	11.4462
Improved Variogram	195.9085	0.7501	7.0059

, and the results have been greatly improved.

5. Experiments

To verify the effectiveness of this method, the aeromagnetic survey data of the Canadian Geological Survey (data sources: Kluane Lake West Aeromagnetic Survey, Residual Total Magnetic Field, NTS 115G/12 and parts of NTS 115G/11, 13, 14 and NTS 115F/9 and 16, Yukon) [18,19] are used for analysis and research. There are abundant deposits and significant magnetic anomalies in the Kluane area. The regional resolution of the geomagnetic map is 50 × 50 m, as shown in Figure 12a.

The geomagnetic data with the distance of 500 m in the line-spacing long the Y direction is taken from the geomagnetic map in Figure 12a, and superimposes the Gaussian white noise with a mean of 0 nT and a variance of 0.1, which is taken as the aeromagnetic acquisition data to be interpolated. The spherical model and the method proposed in this paper are respectively used to carry out Kriging interpolation composition. The interpolation results and error distribution are shown in Figure 13b,e. As can be seen from Figure 12, the trend of the geomagnetic map is about 10° downward along the X direction. Through the analysis of the direction–range diagram, the coordinate azimuth angle 162° is selected in this paper for geometric anisotropy processing, and the comparison results between the processed geomagnetic map and the interpolation map are shown in Figure 13e. The 11 lines in Figure 13 are increased according to the algorithm proposed in this paper (as shown in Figure 13d). It can be seen from the comparison chart that the proposed algorithm is better. Meanwhile, compared with the maximum absolute error (ME), mean absolute error (MAE) and root-mean-square error (RMSE) (as shown in

Table 5. Comparison of Kriging interpolation errors for the aeromagnetic survey data using different variogram models.

Variogram Model	ME (nT)	MAE (nT)	RMSE (nT)
Spherical Variogram	131.7929	10.0630	17.4667
Improved Variogram	92.8156	7.7730	14.4882

), the calculated results are all improved to a certain extent compared with the results without geometric anisotropy.

6. Conclusions

The main work of this paper includes:

• The influence of geometric anisotropy caused by the sparsity of geomagnetic data on the variability of geomagnetic map construction is analyzed;
• The geometric anisotropy correction formula of the variogram of sparse data is derived;
• A geomagnetic interpolation algorithm based on geometric anisotropy correction of variogram is proposed.

At the same time, through the simulation and measured data test, it is verified that the interpolation accuracy of this method has obvious advantages, which greatly reduces the maximum absolute error, average absolute error and root mean square error, and the performance of ME, MAE and RMSE is improved by at least 10%. Furthermore, the performance is improved by four times in the simulation scene with uncomplicated geomagnetic field based on two magnetic dipoles. To sum up, this modified interpolation algorithm can be applied to the construction of local geomagnetic anomaly maps with high precision and high resolution.