# Application of Kriging Algorithm Based on ACFPSO in Geomagnetic Data Interpolation

Authors: Zhijian Zhou; Meng Zhang; Yanzhang Wang; Chao Wang; Ming Ma
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

High-precision geomagnetic field model is the key to magnetic anomaly detection and localization technology. The model is usually constructed through Kriging interpolation. Aiming at the problem of insufficient fitting of variogram in the existing Kriging interpolation methods, this paper proposes a particle swarm optimization algorithm with an adaptive compression factor (ACFPSO). The algorithm utilizes the degree of particle aggregation and the number of iterations to dynamically change the compression factor so as to achieve an effective balance between global optimization and local exploration. The cross-validation results show that the ACFPSO algorithm has the same convergence speed as the conventional particle swarm optimization algorithm, but the convergence accuracy is higher. Compared with the commonly used high-efficiency interpolation methods, such as the plain Kriging, the inverse distance weighting, and the radial basis function, the ACFPSO-optimized Kriging method achieves better performance (the mean absolute error is around 0.3 nT).


---

## Body

## 1. Introduction

Geomagnetic data provide important information for near-surface detection [1]. However, existing measurement methods cannot realize fine-grained measurements on the geomagnetic field. The spatial interpolation algorithm can effectively increase the density of measurements [2]. When traditional interpolation methods are applied to continuous geomagnetic field modelling, limited or sparse geomagnetic data can lead to bias and inconsistent spatial estimation in spatial prediction results. Kriging interpolation as a geostatistical method is one of the widely used interpolation methods in the field of precipitation prediction [3, 4], soil property estimation [5, 6], groundwater statistics [7], and digital terrain modelling [8]. Kriging interpolation may overcome the limitations encountered with current interpolation methods and enhance the quality of geomagnetic data.Fitting the variogram model as a key link in Kriging interpolation usually determines the parameters of the variogram by selecting and fitting the theoretical model of the variogram [9]. Many scholars have made different improvements to the fitting methods in order to improve the accuracy of interpolation. However, the mainstream methods of fitting variograms have their limitations. For example, the least-squares method [10] often falls into local optimum, and the genetic algorithm [11] is of a complicated structure and has more parameters to tune, even though its accuracy is relatively high.In recent years, with the development of nonlinear models, it is not necessary to convert it to a linear model when solving the nonlinear model [12]. The particle swarm optimization (PSO) algorithm as a kind of typical intelligent algorithm has the advantage of finding a globally optimal solution [13]. It can not only directly solve the nonlinear equation but also the variation equation [14]. The PSO algorithm has been widely used for nonlinear fitting due to its fast convergence speed, a small number of hyperparameters, and compatibility with nondifferentiable functions [15, 16]. However, the conventional PSO algorithm is prone to fall into local optimum when solving more complicated problems [17–19]. Therefore, this paper proposes a particle swarm optimization algorithm with an adaptive compression factor and applies an optimized algorithm to Kriging interpolation.


## 2. Kriging Interpolation Based on Optimized Particle Swarm Optimization

### 2.1. Kriging Interpolation

Kriging interpolation is a method based on variogram theory and structural analysis to unbiased optimal estimation of regionalized variables in a finite region [20]. Suppose regionalized variable Zx meets the second-order stationarity conditions. There are n reference points around the interpolation point x0 with attribute value xii=1,2,…,n. A Kriging interpolation model is defined as follows [21]:(1)Z′x0=∑i=1nλiZxi,where λi is the weighting factor of the reference point xi, which is determined by the calculation result of the variogram based on the unbiased and minimum variance conditions. Suppose the two points are separated by h and the function values of two points are Zx and Zx+h, the variogram can be expressed as(2)γh=12VarZx−Zx+h,where h is the lag distance and γh is the experimental variogram.In applied geostatistics, experimental variograms are approximated by the theoretical variogram models [22]. Some commonly used theoretical models mainly include the exponential, the spherical, and the Gaussian models [23]. The model that more often describes the reality geological statistics [24, 25] is the Gaussian model based on the following equation:(3)rh=0,h<0,C0+C1−e−h2/a2,h>0,where C0 is the nugget, a is the range, and C0+C is the sill. Semivariogram modelling and estimation are extremely important for structural analysis and spatial interpolation.Generally, the variogram of a regionalized variable is unknown, and the experimental variogram should be calculated using the data gathered from the regionalized variable [14]. Under the condition that the regionalized variable satisfies the second-order stationary, Kriging equations can be derived from the following formula:(4)∑i=1nλi=1,∑i=1nλiγxi,xj+μ=γxi,x0,where γxi,xj is the variation functional value between the measurement points xi and xj, γxi,x0 is the variation functional value between the observation point xi and the prediction point x0, and μ is the Lagrange multiplier associated with minimizing variance.Finally, the interpolation model can be expressed as follows:(5)γx1,x1…γx1,xn1…………γxn,x1…γxn,xn11……1λ1…λnμ=γx1,x0…γxn,x01.From equation (5), the value of the weight λi is calculated. The estimated value Z′x0 at the predicted point x0 can be obtained by taking formula (1).


### 2.2. Particle Swarm Optimization

PSO is an optimization algorithm based on a population-based global search. It corrects individual mobility strategies through information sharing between populations and individual experience summarization and finally obtains solutions to optimization problems [26]. In the PSO algorithm, suppose the size of the population is N, and the dimension of the solution space is D-dimensional. Xi=Xi1,Xi2,…,XiD is the current position of the particle. Vi=Vi1,Vi2,…,ViD is the current speed of the particle. The fitness value of each particle is denoted as fXi, and the position of the particle is judged according to the fitness value of the particle. The optimal position p best experienced by the particle itself is pbi=pbi1,pbi2,…,pbiD. The optimal position g best of the current group experience is gbi=gbi1,gbi2,…,gbiD [27]. Mathematical expressions of update speed and position of the particle are as follows:(6)vidk+1=vidk+c1R1pbidk−xidk+c2R1gbidk−xidk,xidk+1=xidk+vidk+1,where R1andR2 are the random numbers between 0 and 1, k is the number of iterations, and c1 and c2 are the acceleration coefficients.In order to control and constrain the flight speed of particles, Clerc proposed the concept of compression factor particle swarm optimization (CFPSO) in 1999 [28]. The speed update formula is as follows:(7)vik+1=φvik+c1R1pbidk−xidk+c2R2gbidk−xidk,(8)φ=22−ρ−ρ2−ρ,ρ=ρ1+ρ2,ρ>4,where φ is the compression factor. The compression factor in equation (7) remains unchanged during the algorithm optimization process, but the size of the compression factor directly affects the convergence of the algorithm. The compression factor is too large, and the accuracy of solving the global optimal solution will be reduced. On the contrary, the diversity of the population will be reduced, and the precocity phenomenon will occur [29].


### 2.3. Modified PSO

In order to improve the precision of the CFPSO algorithm, this paper proposes an ACFPSO algorithm, which adaptively adjusts the compression factor according to the search progress of the algorithm, thus improving the ability of the algorithm to search for the global optimal solution.In neural networks, the most commonly used constructive neuron activation function is the sigmoid function [30] as shown in the following equation:(9)Sx=11+ex.This function shows a better balance between linear and nonlinear behaviors [31]. According to the strategy of constructing nonlinear decreasing inertia weight in the literature [32], combined with the smooth transition between the linear and nonlinear functions of the neuron activation function sigmoid, this paper proposes a sigmoid-based compression factor:(10)φ=11+eT/T−ti−c,where T is the maximum number of iterations, ti is the current number of iterations, and c is a constant. The compression factor is large at the beginning, which makes the algorithm have a certain exploration ability. As the number of iterations increases, the compression factor decreases gradually, which makes the algorithm to have local contraction ability, thus improving the ability to search for the global optimal solution.Aggregation also affects the results of optimization in the process of optimization. Therefore, this paper proposes to combine the influence of the number of iterations with the degree of aggregation to realize the dynamic adjustment of the compression factor. The optimization scheme proposed in this paper utilizes average aggregation distance and maximum aggregation distance [33], which formula is as follows:(11)Mean=∑i=1m∑d=1DPid−xid2m,Max=maxi=1,2,…,m∑i=1DPid−xid2,where Mean is the average gathering distance, Max is the maximum gathering distance, m is the total number of particles, D is the dimension of particles, pbi is the optimal position currently searched by particle swarm, and xid is the optimal position currently searched by each particle.When Mean is large, it indicates that the particles are dispersed, and it is necessary to increase the compression factor, thereby enhancing the global search ability of the particles. When Mean is decreasing, it indicates that the particles are converging, and it is necessary to reduce the compression factor, thereby enhancing the local search ability of the particles. In order to get a better grasp of the optimization process of particles, this paper defines the rate of change of particle aggregation distance, which formula is as follows:(12)Mi=MeanMax.The average aggregation distance and the maximum aggregation distance can be obtained once per iteration, and the rate of change of the aggregation distance of this iteration is obtained by calculation.In order to better balance the speed and accuracy of the PSO algorithm, we use the average valueM of the particle aggregation change rate at half the iteration as the criterion, as shown in the following equation:(13)M=∑i=1T/2MiT/2.Whenti<T/2, Mi is too large, indicating that the initial particles are relatively dispersed, and it is necessary to increase φ to make the particles bias toward global convergence and speed up the search. When ti<T/2, it is necessary to reduce φ to bias the particles towards local search. However, when ti<T/2 and Mi>M, the magnitude of w reduction should be increased to make the particles reach local convergence quickly. According to the above analysis, the speed update formula of the ACFPSO algorithm proposed in this paper is as follows:(14)vik+1=φvik+c1R1pbidk−xidk+c2R2gbidk−xidk,(15)φ=k1×11+eT/T−ti−c,ti<T2,k2×11+eT/T−ti−c,ti>T2,Mi<M,k3×11+eT/T−ti−c,ti>T2,Mi>M,where φ is the compression factor, T is the maximum number of iterations, ti is the current number of iterations, and c is a constant. k1, k2, and k3 are the control factors.The compression factor of the entire optimization process changes with the number of iterations and the rate of change of the aggregation distance. When the degree of aggregation of particles is low during the iterative process, the compression factor should be increased so that the poor particles could adaptively increase the search space and prevent the aggregation from rising excessively. When the degree of particle aggregation is high, the compression factor should be reduced, ensuring that the optimal particle performs an accurate search for the surrounding neighborhood. The flow chart of the proposed algorithm is shown in Figure1.Figure 1
The flow chart of the ACFPSO algorithm.


## 2.1. Kriging Interpolation

Kriging interpolation is a method based on variogram theory and structural analysis to unbiased optimal estimation of regionalized variables in a finite region [20]. Suppose regionalized variable Zx meets the second-order stationarity conditions. There are n reference points around the interpolation point x0 with attribute value xii=1,2,…,n. A Kriging interpolation model is defined as follows [21]:(1)Z′x0=∑i=1nλiZxi,where λi is the weighting factor of the reference point xi, which is determined by the calculation result of the variogram based on the unbiased and minimum variance conditions. Suppose the two points are separated by h and the function values of two points are Zx and Zx+h, the variogram can be expressed as(2)γh=12VarZx−Zx+h,where h is the lag distance and γh is the experimental variogram.In applied geostatistics, experimental variograms are approximated by the theoretical variogram models [22]. Some commonly used theoretical models mainly include the exponential, the spherical, and the Gaussian models [23]. The model that more often describes the reality geological statistics [24, 25] is the Gaussian model based on the following equation:(3)rh=0,h<0,C0+C1−e−h2/a2,h>0,where C0 is the nugget, a is the range, and C0+C is the sill. Semivariogram modelling and estimation are extremely important for structural analysis and spatial interpolation.Generally, the variogram of a regionalized variable is unknown, and the experimental variogram should be calculated using the data gathered from the regionalized variable [14]. Under the condition that the regionalized variable satisfies the second-order stationary, Kriging equations can be derived from the following formula:(4)∑i=1nλi=1,∑i=1nλiγxi,xj+μ=γxi,x0,where γxi,xj is the variation functional value between the measurement points xi and xj, γxi,x0 is the variation functional value between the observation point xi and the prediction point x0, and μ is the Lagrange multiplier associated with minimizing variance.Finally, the interpolation model can be expressed as follows:(5)γx1,x1…γx1,xn1…………γxn,x1…γxn,xn11……1λ1…λnμ=γx1,x0…γxn,x01.From equation (5), the value of the weight λi is calculated. The estimated value Z′x0 at the predicted point x0 can be obtained by taking formula (1).


## 2.2. Particle Swarm Optimization

PSO is an optimization algorithm based on a population-based global search. It corrects individual mobility strategies through information sharing between populations and individual experience summarization and finally obtains solutions to optimization problems [26]. In the PSO algorithm, suppose the size of the population is N, and the dimension of the solution space is D-dimensional. Xi=Xi1,Xi2,…,XiD is the current position of the particle. Vi=Vi1,Vi2,…,ViD is the current speed of the particle. The fitness value of each particle is denoted as fXi, and the position of the particle is judged according to the fitness value of the particle. The optimal position p best experienced by the particle itself is pbi=pbi1,pbi2,…,pbiD. The optimal position g best of the current group experience is gbi=gbi1,gbi2,…,gbiD [27]. Mathematical expressions of update speed and position of the particle are as follows:(6)vidk+1=vidk+c1R1pbidk−xidk+c2R1gbidk−xidk,xidk+1=xidk+vidk+1,where R1andR2 are the random numbers between 0 and 1, k is the number of iterations, and c1 and c2 are the acceleration coefficients.In order to control and constrain the flight speed of particles, Clerc proposed the concept of compression factor particle swarm optimization (CFPSO) in 1999 [28]. The speed update formula is as follows:(7)vik+1=φvik+c1R1pbidk−xidk+c2R2gbidk−xidk,(8)φ=22−ρ−ρ2−ρ,ρ=ρ1+ρ2,ρ>4,where φ is the compression factor. The compression factor in equation (7) remains unchanged during the algorithm optimization process, but the size of the compression factor directly affects the convergence of the algorithm. The compression factor is too large, and the accuracy of solving the global optimal solution will be reduced. On the contrary, the diversity of the population will be reduced, and the precocity phenomenon will occur [29].


## 2.3. Modified PSO

In order to improve the precision of the CFPSO algorithm, this paper proposes an ACFPSO algorithm, which adaptively adjusts the compression factor according to the search progress of the algorithm, thus improving the ability of the algorithm to search for the global optimal solution.In neural networks, the most commonly used constructive neuron activation function is the sigmoid function [30] as shown in the following equation:(9)Sx=11+ex.This function shows a better balance between linear and nonlinear behaviors [31]. According to the strategy of constructing nonlinear decreasing inertia weight in the literature [32], combined with the smooth transition between the linear and nonlinear functions of the neuron activation function sigmoid, this paper proposes a sigmoid-based compression factor:(10)φ=11+eT/T−ti−c,where T is the maximum number of iterations, ti is the current number of iterations, and c is a constant. The compression factor is large at the beginning, which makes the algorithm have a certain exploration ability. As the number of iterations increases, the compression factor decreases gradually, which makes the algorithm to have local contraction ability, thus improving the ability to search for the global optimal solution.Aggregation also affects the results of optimization in the process of optimization. Therefore, this paper proposes to combine the influence of the number of iterations with the degree of aggregation to realize the dynamic adjustment of the compression factor. The optimization scheme proposed in this paper utilizes average aggregation distance and maximum aggregation distance [33], which formula is as follows:(11)Mean=∑i=1m∑d=1DPid−xid2m,Max=maxi=1,2,…,m∑i=1DPid−xid2,where Mean is the average gathering distance, Max is the maximum gathering distance, m is the total number of particles, D is the dimension of particles, pbi is the optimal position currently searched by particle swarm, and xid is the optimal position currently searched by each particle.When Mean is large, it indicates that the particles are dispersed, and it is necessary to increase the compression factor, thereby enhancing the global search ability of the particles. When Mean is decreasing, it indicates that the particles are converging, and it is necessary to reduce the compression factor, thereby enhancing the local search ability of the particles. In order to get a better grasp of the optimization process of particles, this paper defines the rate of change of particle aggregation distance, which formula is as follows:(12)Mi=MeanMax.The average aggregation distance and the maximum aggregation distance can be obtained once per iteration, and the rate of change of the aggregation distance of this iteration is obtained by calculation.In order to better balance the speed and accuracy of the PSO algorithm, we use the average valueM of the particle aggregation change rate at half the iteration as the criterion, as shown in the following equation:(13)M=∑i=1T/2MiT/2.Whenti<T/2, Mi is too large, indicating that the initial particles are relatively dispersed, and it is necessary to increase φ to make the particles bias toward global convergence and speed up the search. When ti<T/2, it is necessary to reduce φ to bias the particles towards local search. However, when ti<T/2 and Mi>M, the magnitude of w reduction should be increased to make the particles reach local convergence quickly. According to the above analysis, the speed update formula of the ACFPSO algorithm proposed in this paper is as follows:(14)vik+1=φvik+c1R1pbidk−xidk+c2R2gbidk−xidk,(15)φ=k1×11+eT/T−ti−c,ti<T2,k2×11+eT/T−ti−c,ti>T2,Mi<M,k3×11+eT/T−ti−c,ti>T2,Mi>M,where φ is the compression factor, T is the maximum number of iterations, ti is the current number of iterations, and c is a constant. k1, k2, and k3 are the control factors.The compression factor of the entire optimization process changes with the number of iterations and the rate of change of the aggregation distance. When the degree of aggregation of particles is low during the iterative process, the compression factor should be increased so that the poor particles could adaptively increase the search space and prevent the aggregation from rising excessively. When the degree of particle aggregation is high, the compression factor should be reduced, ensuring that the optimal particle performs an accurate search for the surrounding neighborhood. The flow chart of the proposed algorithm is shown in Figure1.Figure 1
The flow chart of the ACFPSO algorithm.


## 3. Discussion

In order to test the performance of the ACFPSO algorithm, this paper addresses three benchmarking problems (as shown in Table1) for comparison.Table 1
Test functions and their properties.
Test function
Formula
Ranges
Theoretical optimal
Sphere
f
x
=
∑
i
=
1
n
x
i
2
(−100, 100)
0
Rastrigin
f
x
=
∑
i
=
1
n
x
i
2
−
10
cos
2
π
x
i
+
10
(−100, 100)
0
Griewangk
f
x
=
1
4000
∑
i
=
1
n
x
i
2
−
∏
i
=
1
n
cos
x
i
/
i
+
1
(−100, 100)
0The above three test problems are used to compare the performance of PSO, CFPSO, and ACFPSO algorithms. In the experiment, the number of particles is set to 40, the number of dimensions is set to 5 and 10, the number of iterations is set to 100 and 500, and the acceleration coefficientsc1=c2=2.05 [17]. In equation (15), c is a constant, and k1, k2 and k3 are the control factors. Considering the convergence accuracy and speed of the ACFPSO algorithm, the good optimization results are obtained when c=2.5, k1=1.5, k2=0.8, and k3=0.2. By setting different dimensions and the number of iterations to verify the performance of the algorithm, each algorithm runs independently for 30 times, and the average value, standard deviation, and average time consumption of the function optimal fitness value are used as evaluation indexes. The computing unit used for the test is an Intel Core i5-6400 CPU (2.70 GHz), and the soft environment is Octave. The test results are shown in Tables 2–5.Table 2
Test results with 5 dimensions and 100 iterations.
Test function
Index
PSO
CFPSO
ACFPSO
Sphere
Mean
0.0073
2.0150e − 04
7.7923e − 09
Deviation
0.0062
2.8657e − 04
2.2586e − 08
Time (s)
0.027223
0.022201
0.028824
Rastrigin
Mean
0.0112
2.1919e − 04
2.3561e − 08
Deviation
0.0081
4.7060e − 04
4.7727e − 08
Time (s)
0.028003
0.025485
0.026312
Griewangk
Mean
0.0085
9.3705e − 05
1.4708e − 09
Deviation
0.0063
1.3605e − 04
2.2445e − 09
Time (s)
0.027149
0.020732
0.024963Table 3
Test results with 10 dimensions and 100 iterations.
Test function
Index
PSO
CFPSO
ACFPSO
Sphere
Mean
0.1516
3.2353e − 04
3.3724e − 07
Deviation
0.0727
4.4511e − 04
7.7252e − 07
Time (s)
0.022765
0.020077
0.028838
Rastrigin
Mean
0.1706
2.3591e − 04
3.6331e − 07
Deviation
0.0748
3.2564e − 04
7.2917e − 07
Time (s)
0.018033
0.020819
0.026974
Griewangk
Mean
0.0200
6.2005e − 05
2.5018e − 07
Deviation
0.0084
1.2158e − 04
6.2321e − 07
Time (s)
0.020094
0.024233
0.024334Table 4
Test results with 5 dimensions and 500 iterations.
Test function
Index
PSO
CFPSO
ACFPSO
Sphere
Mean
0.0066
2.5869e − 24
0
Deviation
0.0047
1.4169e − 23
0
Time (s)
0.112825
0.091152
0.108400
Rastrigin
Mean
0.0059
2.0724e − 15
0
Deviation
0.0038
8.3850e − 15
0
Time (s)
0.090972
0.083736
0.088310
Griewangk
Mean
0.0012
9.6219e − 17
0
Deviation
8.3142e − 04
1.9059e − 16
0
Time (s)
0.090130
0.092149
0.096418Table 5
Test results with 10 dimensions and 500 iterations.
Test function
Index
PSO
CFPSO
ACFPSO
Sphere
Mean
0.00375
2.39e − 25
1.53e − 35
Deviation
0.00259
1.31e − 24
4.46e − 35
Time (s)
0.2457
0.1154
0.1175
Rastrigin
Mean
0.00225
1.77e − 16
0
Deviation
0.00152
7.15e − 16
0
Time (s)
0.2375
0.1191
0.1113
Griewangk
Mean
0.00928
7.40e − 18
0
Deviation
0.00429
4.05e − 17
0
Time (s)
0.2391
0.1157
0.1183It can be seen from Tables2–5 that under the same parameter settings, the ACFPSO algorithm has obvious advantages over the PSO and the CFPSO algorithms in terms of the mean and standard deviation of optimal fitness values for the three benchmark problems. The optimization speeds of CFPSO and ACFPSO algorithms are basically the same, but both speeds are much higher than those of the PSO algorithm.The iterative evolution curves of three different particle swarm optimization algorithms are compared in Figures2–5. It can be seen that the ACFPSO algorithm and the CFPSO algorithm have greatly expedited convergence speed and improved convergence accuracy compared with the PSO algorithm. ACFPSO and CFPSO have the same convergence speed (20 iterations). However, the comparison between Tables 2–5 shows that ACFPSO has higher accuracy than CFPSO. Therefore, the ACFPSO algorithm proposed in this paper has better dynamic optimization performance and can be used as an improvement of the PSO algorithm.Figure 2
Evolution curves with 5 dimensions and 100 iterations. (a) Sphere function. (b) Rastrigin function. (c) Griewangk function.
(a)
(b)
(c)Figure 3
Evolution curves with 10 dimensions and 100 iterations. (a) Sphere function. (b) Rastrigin function. (c) Griewangk function.
(a)
(b)
(c)Figure 4
Evolution curves with 5 dimensions and 500 iterations. (a) Sphere function. (b) Rastrigin function. (c) Griewangk function.
(a)
(b)
(c)Figure 5
Evolution curves with 10 dimensions and 500 iterations. (a) Sphere function. (b) Rastrigin function. (c) Griewangk function.
(a)
(b)
(c)


## 4. Experiments on Real Data

Section3 has validated that the ACFPSO algorithm has higher convergence precision and can be used as a method to fit the Kriging variogram. In order to test the accuracy and reliability of the ACFPSO-Kriging interpolation algorithm, this paper uses the measured data to verify the optimized interpolation algorithm.


### 4.1. Data Acquisition

The data for the test are a set of aeromagnetic data, which were acquired by using the measurement method of the helicopter towed bird, as shown in Figure6(a). This towed bird is connected to the helicopter by a streamer with a length of 50 m. The flight altitude ranges from 100 m to 200 m. A total of 25 survey lines cover the area as shown in Figure 6(b). The length of each survey line is around 3 km, and the line spacing is around 100 m. 20, 000 samples of the geomagnetic intensity were acquired during the flight at a sampling rate of 10 Hz. Of these data, 10 percent is taken out as the reference, and the other 90 percent is to be interpolated. The interpolated samples are compared with the ground truth to evaluate the proposed method’s performance.Figure 6
Aeromagnetic data acquisition diagram and data acquisition path.
(a)
(b)


### 4.2. Comparison of Fitting Performance

The geostatistical analysis uses the semivariogram to quantify the spatial variation of a regionalized variable and derive important parameters used for Kriging spatial interpolation. Before fitting the variogram, it is necessary to verify the effectiveness of the proposed algorithm for adaptively adjusting the compression factor based on the number of iterations and the rate of change of the aggregation distance. The parameter calculation of the variogram model is performed using the geomagnetic data measured above, and the variogram model is a Gaussian model. In the experiment, the number of particles is 40, the number of dimensions is 3, the number of iterations is 500, and the acceleration coefficientsc1=c2=2.05. Figures 7(a) and 7(b) are comparisons of particle aggregation distance rates between CFPSO and ACFPSO algorithms for solving variogram model parameters.Figure 7
Comparison of the change rate of aggregation distance. (a) Based on the CFPSO algorithm. (b) Based on the ACFPSO algorithm.
(a)
(b)It can be seen that the ACFPSO algorithm proposed in this paper is more accurate and faster in solving the parameters of the variogram model. Then the parameters of the variogram model are solved by cross-validation. The relevant parameters of the Gaussian model were obtained by fitting the experimental variograms by LS, PSO, CFPSO, and ACFPSO algorithms, as shown in Table6. The comparison curves of the fitted curves are shown in Figure 8.Table 6
Theoretical variogram model parameter value.
a
C
C
0
LS
3430.3482
19696.2216
1010.3708
PSO
4281.4297
15063.7659
1124.0739
CFPSO
6714.7519
16020.3347
1017.7951
ACFPSO
2323.76524
16038.0072
933.5771Figure 8
Comparison between ACFPSO and other three algorithms in the fitting variogram.
(a)
(b)
(c)It can be seen from Figure8 that the fitting effect by the ACFPSO algorithm is closer to the trend of the experimental variogram value than the least-squares method and the PSO and the CFPSO algorithms. Especially when the lag distance is more than 2000 m, the ACFPSO algorithm fits more closely with the experimental variogram curve. In this paper, the mean absolute error (MAE) and the mean absolute percentage error (MAPE) are introduced to evaluate the four fitting methods. The calculation formulas are(16)MAE=1N∑I=1Nzi−zi∗,MAPE=1N∑I=1Nzi−zi∗zi.Figure9 is a histogram of the distribution of MAE and MAPE. It can be seen from Figure 9 that the ACFPSO algorithm has the highest accuracy in fitting the variogram.Figure 9
Distribution histograms of (a) MAE and (b) MAPE.
(a)
(b)


### 4.3. Comparison and Analysis of Interpolation Effects

In order to verify the better characteristics of the optimized Kriging interpolation algorithm, the more efficient interpolation methods such as ordinary Kriging, the inverse distance weighting (IDW), and the radial basis function (RBF) are compared with the interpolation results of the optimized Kriging interpolation algorithm proposed in this paper. The evaluation parameters were the mean error (ME), the mean absolute error (MAE), the root-mean-square error (RMSE), andR square (R2). The formulas for evaluating the parameters are as follows:(17)ME=1N∑i=1Nzi−zi∗,MSE=1N∑i=1Nzi−zi∗2,RMSE=1N∑i=1Nzi−zi∗2,R2=1−∑i=1NZi−Zi∗2∑i=1Nzi−z¯i2,where zi is the measured true value, zi∗ is the predicted value, z¯i is the measured mean value, and N is the verification point number. The error analysis is performed on the interpolation results using ME, MAE, MSE, RMSE, R2, and time (s), as shown in Table 7.Table 7
Interpolation error analysis of different methods.
ME
MAE
RMSE
R
2
Time (s)
IDW
0.1053
2.2381
3.4090
0.9984
3.4934
RBF
0.0861
0.9171
1.8177
0.9995
5.4882
LS-Kriging
0.0823
0.5157
1.0966
0.9998
12.2869
PSO-Kriging
0.0726
0.4067
0.9743
0.9999
12.0434
CFPSO-Kriging
0.0716
0.3883
0.9534
0.9999
11.9235
ACFPSO-Kriging
0.0706
0.2949
0.8487
0.9999
11.9621It can be seen from Table7 that four different Kriging interpolations have better accuracy compared with the IDW and the RBF methods. Kriging interpolation based on the ACFPSO algorithm outperforms the LS-Kriging, the PSO-Kriging, and the CFPSO-Kriging.The above method is used to obtain the theoretical variogram model for Kriging interpolation of predicted points. Figures10(a) and 10(b) are, respectively, contour maps of the geomagnetic field measured data and the geomagnetic field data formed by Kriging interpolation based on ACFPSO, and the white squares in the two figures are magnetic anomalies. It can be seen from the comparison between Figures 10(a) and 10(b) that the trend of the measured magnetic field data and several magnetic anomalies can be accurately reflected. Therefore, Kriging interpolation based on the ACFPSO algorithm is effective in the field of magnetic anomaly detection.Figure 10
Geomagnetic data contour map. (a) Actual measurement data. (b) Interpolated data based on the ACFPSO’s Kriging method.
(a)
(b)


## 4.1. Data Acquisition

The data for the test are a set of aeromagnetic data, which were acquired by using the measurement method of the helicopter towed bird, as shown in Figure6(a). This towed bird is connected to the helicopter by a streamer with a length of 50 m. The flight altitude ranges from 100 m to 200 m. A total of 25 survey lines cover the area as shown in Figure 6(b). The length of each survey line is around 3 km, and the line spacing is around 100 m. 20, 000 samples of the geomagnetic intensity were acquired during the flight at a sampling rate of 10 Hz. Of these data, 10 percent is taken out as the reference, and the other 90 percent is to be interpolated. The interpolated samples are compared with the ground truth to evaluate the proposed method’s performance.Figure 6
Aeromagnetic data acquisition diagram and data acquisition path.
(a)
(b)


## 4.2. Comparison of Fitting Performance

The geostatistical analysis uses the semivariogram to quantify the spatial variation of a regionalized variable and derive important parameters used for Kriging spatial interpolation. Before fitting the variogram, it is necessary to verify the effectiveness of the proposed algorithm for adaptively adjusting the compression factor based on the number of iterations and the rate of change of the aggregation distance. The parameter calculation of the variogram model is performed using the geomagnetic data measured above, and the variogram model is a Gaussian model. In the experiment, the number of particles is 40, the number of dimensions is 3, the number of iterations is 500, and the acceleration coefficientsc1=c2=2.05. Figures 7(a) and 7(b) are comparisons of particle aggregation distance rates between CFPSO and ACFPSO algorithms for solving variogram model parameters.Figure 7
Comparison of the change rate of aggregation distance. (a) Based on the CFPSO algorithm. (b) Based on the ACFPSO algorithm.
(a)
(b)It can be seen that the ACFPSO algorithm proposed in this paper is more accurate and faster in solving the parameters of the variogram model. Then the parameters of the variogram model are solved by cross-validation. The relevant parameters of the Gaussian model were obtained by fitting the experimental variograms by LS, PSO, CFPSO, and ACFPSO algorithms, as shown in Table6. The comparison curves of the fitted curves are shown in Figure 8.Table 6
Theoretical variogram model parameter value.
a
C
C
0
LS
3430.3482
19696.2216
1010.3708
PSO
4281.4297
15063.7659
1124.0739
CFPSO
6714.7519
16020.3347
1017.7951
ACFPSO
2323.76524
16038.0072
933.5771Figure 8
Comparison between ACFPSO and other three algorithms in the fitting variogram.
(a)
(b)
(c)It can be seen from Figure8 that the fitting effect by the ACFPSO algorithm is closer to the trend of the experimental variogram value than the least-squares method and the PSO and the CFPSO algorithms. Especially when the lag distance is more than 2000 m, the ACFPSO algorithm fits more closely with the experimental variogram curve. In this paper, the mean absolute error (MAE) and the mean absolute percentage error (MAPE) are introduced to evaluate the four fitting methods. The calculation formulas are(16)MAE=1N∑I=1Nzi−zi∗,MAPE=1N∑I=1Nzi−zi∗zi.Figure9 is a histogram of the distribution of MAE and MAPE. It can be seen from Figure 9 that the ACFPSO algorithm has the highest accuracy in fitting the variogram.Figure 9
Distribution histograms of (a) MAE and (b) MAPE.
(a)
(b)


## 4.3. Comparison and Analysis of Interpolation Effects

In order to verify the better characteristics of the optimized Kriging interpolation algorithm, the more efficient interpolation methods such as ordinary Kriging, the inverse distance weighting (IDW), and the radial basis function (RBF) are compared with the interpolation results of the optimized Kriging interpolation algorithm proposed in this paper. The evaluation parameters were the mean error (ME), the mean absolute error (MAE), the root-mean-square error (RMSE), andR square (R2). The formulas for evaluating the parameters are as follows:(17)ME=1N∑i=1Nzi−zi∗,MSE=1N∑i=1Nzi−zi∗2,RMSE=1N∑i=1Nzi−zi∗2,R2=1−∑i=1NZi−Zi∗2∑i=1Nzi−z¯i2,where zi is the measured true value, zi∗ is the predicted value, z¯i is the measured mean value, and N is the verification point number. The error analysis is performed on the interpolation results using ME, MAE, MSE, RMSE, R2, and time (s), as shown in Table 7.Table 7
Interpolation error analysis of different methods.
ME
MAE
RMSE
R
2
Time (s)
IDW
0.1053
2.2381
3.4090
0.9984
3.4934
RBF
0.0861
0.9171
1.8177
0.9995
5.4882
LS-Kriging
0.0823
0.5157
1.0966
0.9998
12.2869
PSO-Kriging
0.0726
0.4067
0.9743
0.9999
12.0434
CFPSO-Kriging
0.0716
0.3883
0.9534
0.9999
11.9235
ACFPSO-Kriging
0.0706
0.2949
0.8487
0.9999
11.9621It can be seen from Table7 that four different Kriging interpolations have better accuracy compared with the IDW and the RBF methods. Kriging interpolation based on the ACFPSO algorithm outperforms the LS-Kriging, the PSO-Kriging, and the CFPSO-Kriging.The above method is used to obtain the theoretical variogram model for Kriging interpolation of predicted points. Figures10(a) and 10(b) are, respectively, contour maps of the geomagnetic field measured data and the geomagnetic field data formed by Kriging interpolation based on ACFPSO, and the white squares in the two figures are magnetic anomalies. It can be seen from the comparison between Figures 10(a) and 10(b) that the trend of the measured magnetic field data and several magnetic anomalies can be accurately reflected. Therefore, Kriging interpolation based on the ACFPSO algorithm is effective in the field of magnetic anomaly detection.Figure 10
Geomagnetic data contour map. (a) Actual measurement data. (b) Interpolated data based on the ACFPSO’s Kriging method.
(a)
(b)


## 5. Conclusions

In this paper, the ACFPSO algorithm is proposed to fit the semivariance function to solve the problem that the conventional Kriging interpolation method does not fit the semivariance function well. The experimental results show that the ACFPSO algorithm has the same convergence speed as the conventional particle swarm optimization algorithm, but the convergence accuracy is higher. Kriging interpolation based on the ACFPSO algorithm outperforms the IDW, the RBF, and the other Kriging interpolation (the MAE is around 0.3 nT). Compared with the conventional Kriging algorithm, the MAE is decreased by 42.82%, the MSE is decreased by 40.10%, and the RMSE is decreased by 22.61%. It can be concluded that the optimized Kriging interpolation algorithm proposed in this paper is feasible and effective in the application of geomagnetic data interpolation.



---

*Source: 1574918-2019-12-21.xml*