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).

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 [1719]. 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 meets the second-order stationarity conditions. There are reference points around the interpolation point with attribute value . A Kriging interpolation model is defined as follows [21]:where is the weighting factor of the reference point , 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 and the function values of two points are and , the variogram can be expressed aswhere is the lag distance and 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:where is the nugget, a is the range, and 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:where is the variation functional value between the measurement points and , is the variation functional value between the observation point and the prediction point , and is the Lagrange multiplier associated with minimizing variance.

Finally, the interpolation model can be expressed as follows:

From equation (5), the value of the weight is calculated. The estimated value at the predicted point 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 , and the dimension of the solution space is D-dimensional. is the current position of the particle. is the current speed of the particle. The fitness value of each particle is denoted as , 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 . The optimal position g best of the current group experience is [27]. Mathematical expressions of update speed and position of the particle are as follows:where are the random numbers between 0 and 1, k is the number of iterations, and and 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: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:

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:where is the maximum number of iterations, is the current number of iterations, and 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: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, is the optimal position currently searched by particle swarm, and 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:

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 value of the particle aggregation change rate at half the iteration as the criterion, as shown in the following equation:

When , 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 , it is necessary to reduce to bias the particles towards local search. However, when and , the magnitude of 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:where is the compression factor, is the maximum number of iterations, is the current number of iterations, and is a constant. , , and 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 Figure 1.

3. Discussion

In order to test the performance of the ACFPSO algorithm, this paper addresses three benchmarking problems (as shown in Table 1) for comparison.

The 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 coefficients [17]. In equation (15), is a constant, and , and are the control factors. Considering the convergence accuracy and speed of the ACFPSO algorithm, the good optimization results are obtained when , , , and . 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 25.

It can be seen from Tables 25 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 Figures 25. 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 25 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.

4. Experiments on Real Data

Section 3 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 Figure 6(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.

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 coefficients . Figures 7(a) and 7(b) are comparisons of particle aggregation distance rates between CFPSO and ACFPSO algorithms for solving variogram model parameters.

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 Table 6. The comparison curves of the fitted curves are shown in Figure 8.

It can be seen from Figure 8 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

Figure 9 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.

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), and R square (). The formulas for evaluating the parameters are as follows:where is the measured true value, is the predicted value, 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, , and time (s), as shown in Table 7.

It can be seen from Table 7 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. Figures 10(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.

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.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the Key Projects of National Key R&D Programs for Major Natural Disaster Monitoring, Early Warning and Prevention (Grant no. 2018YFC1503803).