A Method to Remotely Track a Magnetic Target Using a Scalar Magnetometer Array

The orientation of a vector magnetic sensor can affect the measurement accuracy of magnetic anomaly, thereby increasing the localization error of a magnetic target. Compared with vector magnetic sensor, the measurement of the scalar magnetic sensor is almost not influenced by its orientation. Therefore, we present a method for tracking the magnetic target with a static scalar magnetometer array. In this study, the magnitude of the target’s magnetic moment is a key parameter. We isolate it and formulate an optimization problem based on it to estimate the position and magnetic parameters of the target. To calculate the solution of this optimization problem, a dedicated particle swarm optimization (PSO) algorithm is developed. Then, we define a quality index to evaluate the solution calculated by the optimization problem. The proposed method was validated by the simulation and the real data collected when an SUV car was passing by the array on a straight path. The results show that the tracked trajectory is very close to the true trajectory and the quality index can be used as a criterion to allow accepting or rejecting the localization of the target.


Introduction
A target containing ferromagnetic material can be detected by magnetic sensors due to a magnetic field generated by it. Along with the development of the techniques for the magnetic anomaly detection, researchers have proposed several methods for target localization or tracking [1][2][3][4][5][6][7][8][9]. Usually, we can describe the target with six parameters x, y, z, m x , m y , m z . A vector sensor can provide much valuable target information. In order to calculate six parameters of the target, at least two vector magnetometers are needed to build six nonlinear functions. Thus, a sensor array is widely used to locate the magnetic target. Wynn [1] proposed a method of magnetic dipole localization based on the magnetic gradient. Nara et al. [2] obtained a localization formula for magnetic dipole localization by magnetic vectors and their spatial gradients. In this method, the accuracy of localization is highly sensitive to noise in magnetic anomaly field vector data. The measurement data consist of the ambient geomagnetic field and the magnetic fields produced by the target. Therefore, it is hard to obtain the accurate magnetic anomaly field from vector magnetometers in practice. Generally, the gradient of the ambient geomagnetic field is small, and the gradient tensor of the magnetic measurement can reflect the magnetic anomaly of the target. In order to overcome the effect of the geomagnetic field, magnetic gradient tensor arrays comprising multiple vector magnetometers were proposed to locate the magnetic target [4,5,[10][11][12]. Liu and Wang [10] proposed a method for locating the magnetic target using a 3-axis magnetic sensor array. Gang et al. [5] presented a rotationally invariant method for locating a magnetic target using a magnetic gradient tensor system. Despite the obvious advantages in localization with the gradient tensor, challenges exist with gradient tensor measurements. The reasons are as follows: (1) Rotational vibrations of the array will lead to the orientation error of each vector sensor. (2) The strict alignment for the vector sensors in the array is very difficult. (3) Three axes of the vector sensor are not exactly perpendicular. These reasons will affect the magnetic anomaly measurements [13].
However, a scalar magnetic sensor such as an optically pumped magnetometer has an important advantage in that it is relatively insensitive to its orientation. The optically pumped magnetometer does not work in only two orientations when the optical axis is parallel to the ambient field or perpendicular to the ambient field. Therefore, it can work well for a large range of angles. In many location methods, the Euler deconvolution method [14][15][16], which is important for magnetic data interpretation, is used to determine the depth and location of the static target using a moving sensor. In this method, magnetic field grid data (profile or map) are used to determine the position or boundary [17]. Thus, magnetic field data within the region need to be measured.
In this paper, we focus on using a static scalar magnetometer array instead of a stationary array of vector magnetometers to track a moving magnetic target. We present a method for determining the parameters of the moving target. In this method, we isolate the magnitude of the magnetic moment and formulate an optimization problem based on it to estimate the position and magnetic parameters of the target. This can reduce the number of unknown parameters of the optimization problem. Then, we develop a dedicated PSO algorithm to calculate the solution of this optimization problem without the initial parameter values. The PSO algorithm cannot always converge to a global optimal solution in every calculation. Thus, we define a quality index as the criterion for accepting or rejecting the calculated solution in the PSO algorithm.

Localization Theory Based on Scalar Magnetometer Array
The magnetic target can be considered a magnetic dipole when the distance between it and a magnetic sensor is 3 times larger than the size of the dipole's largest dimension [18]. The magnetic field vector B A generated by the dipole is a function of the distance r between the target and the sensor and the magnetic moment m of its own. We build an array with seven scalar magnetometers. The center of the array is located at the origin of a Cartesian coordinate system with the X-axis toward to the geographical north, the Y-axis toward to the geographical east, and Z-axis downward to the ground, as shown in Figure 1. In the coordinate system, the magnetic field generated by the target can be described as follows: where μ 0 = 4π × 10 -7 H/m is the permeability of free space. r is the distance from the target to the sensor. m is the magnitude of the magnetic moment. θ is the rotation angle from the magnetic moment vector toward the X-Y plane, and φ is the rotation angle from the positive X-axis in the X-Y plane. x = x t − x s , y = y t − y s , and z = z t − z s . In practice, the measured magnetic field, B m , is the vector sum of the ambient geomagnetic field B E and the target magnetic field B A [19]: The magnetic field generated by the remote target is far less than the ambient geomagnetic field at the measurement point. According to the approximation, we can rewrite B m in the scalar form as follows: Hence, the magnetic anomaly ΔB measured by a scalar magnetometer can be given as follows: where u denotes the direction vector of B E . I and D denote the inclination and declination of the ambient geomagnetic field, respectively. It is seen from (4) that the magnetic anomaly can be regarded as the projection of B A on B E .
When the magnitude of the magnetic anomaly generated by the target is small, geomagnetic field variation will affect the localization. Generally, the broad characteristics of geomagnetic field variation are consistent over the local region. We can consider that its effect is the same for the measurement of each magnetometer in the array, shown in Figure 1. Therefore, we can eliminate the effect of geomagnetic field variation through the subtraction with two magnetometers and obtain an equation as follows: where i and j are the indices of individual sensors in the array.
There are six parameters in (5), which can be combined into a single six-dimensional vector as p = x t , y t , z t , θ, φ, m . Equation (5) is a high-order nonlinear function, which usually is solved by optimization methods. We define an objective error function as where ΔB m ij is the magnetic field measured by the ith and jth magnetometers. ΔB c ij p is the magnetic field calculated by (5). It can be seen that (7) represents the mean-squared error between the measured and the calculated magnetic fields. Thus, this is a least square error problem. The parameter vector p can be solved by minimizing the error through an appropriate optimization algorithm. This is the normal method for tracking the magnetic target.
Because (5) is a high-order nonlinear function, there are some local optimal solutions in the solution space. It is known that the objective function can affect the performance of the algorithm. If the function is simpler and has fewer parameters, the accuracy of the solution calculated by the optimization algorithm is higher. Thus, we isolate the magnitude of the magnetic moment and rewrite (5) as According to (8), we define a new objective error function as Thus, the new error function has one less parameter than (7). This is good for improving the accuracy of the solution.

Journal of Sensors
Sometimes, the optimization algorithm cannot jump out from the local optima when trapped in it and the algorithm may fail to give the right global solution. Thus, we should define a criterion to evaluate the solution and determine whether to accept or reject the solution according to the criterion.
Once the parameters x t , y t , z t , θ, φ are calculated by the optimization algorithm, the magnitude of the magnetic moment is easily calculated according to (8). Therefore, we can calculate a series of m according to the different magnetic anomalies ΔB ij , which are measured in the array. Theoretically, if the parameters x t , y t , z t , θ, φ of the target are the true value, the calculated moments m are equal to the same value. However, the measurement error leads to the error of parameters x t , y t , z t , θ, φ calculated by the method. Thus, it is highly unlikely for the calculated moments, m, to have the same value. When the parameter errors are small, the calculated moments m are close to the same value. Conversely, when the error of parameters x t , y t , z t , θ, φ is large, the fluctuation of calculated moments m is also large. Therefore, according to the magnetic moment, we define the criterion for the quality of the solution called the quality index as The magnitude of the magnetic moment m with different subscripts (p or q) denotes the calculated magnitude of the magnetic moment using the different ΔB ij .
It is noted that the expression ∑ p,q m p /m q − 1 2 is the squared relative error between the calculated magnitudes of the magnetic moment using the different ΔB ij . The value of this expression is in the range 0, ∞ . It is more convenient for the quality index to be in the interval 0, 1 . Thus, we introduce the exponential function and give the expression of the quality index as (10).
From (10), the closer the values of m calculated by different ΔB are, the closer the expression QI is to 1. This means that the parameters x t , y t , z t , θ, φ of the target calculated by the PSO algorithm are close to the true values. On the contrary, if the parameters x t , y t , z t , θ, φ of the target calculated by the PSO algorithm are the local optimal solution, the fluctuation of calculated moments m is also large. This can cause the value of QI to be far from 1. Therefore, the quality index can be used as an indicator for determining whether the calculated solution is a global optimal solution.

Localization Algorithm Based on the PSO Algorithm
There are many conventional algorithms (such as Levenberg-Marquardt, Newton-Gauss algorithm, and Fletcher-Powell), which can be used to solve the optimization problem. The Levenberg-Marquardt (L-M) algorithm can be regarded as a combination of the steepest descent method and the Newton-Gauss method. The main advantage of the L-M algorithm is its rapidity and its well-known and reliable implementation [20]. However, the L-M algorithm cannot guarantee convergence to a global minimum unless the initial parameter guesses are appropriate [21]. In the tracking system, it is very difficult to give a relatively good initial value of the parameters. The PSO algorithm, which works on the social behavior of particles in a swarm [22], can be used to solve the optimization problem without the initial parameter guesses. We only need to provide the range of solution space (or the range of the parameters). All the particles are considered a swarm, and each particle flies around in the solution space. The particle's position is dynamically updated according to its own experience and the swarm's experience in the iterative process [23,24]. The advantages of the PSO algorithm are that it converges rapidly and it has few adjustable parameters. It should be noted that the optimization algorithm may converge to the local minimum and cannot give the optimal solution. Thus, we define the quality index, which is used as an indicator for determining whether the calculated solution is a global optimal solution. In order to better understand the PSO algorithm, the detailed descriptions of some key terms in the PSO algorithm are given as follows.
(1) Particle flies in the solution space. It has two attributes, position and velocity.
(2) Position denotes the candidate solution for the optimization problem. It is dynamically updated for finding the optimal solution according to the velocity.
(3) Velocity determines the flight direction and distance of the particle in the next step. It is dynamically updated according to the particle's best position and the swarm's best position.
(4) Fitness function is a function for describing optimization problem. It is used to evaluate the position of each particle.
(5) Fitness is the value of the fitness function with the particle's position. It is used to reflect the quality of the position in the optimization.
More details of the PSO algorithm can be seen in [22,25,26]. Here, we develop a dedicated PSO algorithm for tracking the target based on a scalar magnetometer array and give the pseudocode of the algorithm.
In the algorithm, the parameter w called inertia weight is used to bring about a balance between the global and local search. When the value of w is larger, the algorithm is easier to find the global optimal solution, but converges slowly. On the contrary, when the value of w is smaller, the algorithm converges faster, but is easier to fall into local optimal solution area. Thus, the value of w is dynamically adjusted during the iteration. The parameters c 1 and c 2 , called acceleration factors, are used to control the flight of the particle to the particle's best position and the swarm's best position. Researchers have done much work on how to select the ideal parameters for PSO implementation [27][28][29]. In [28,29], inertia weight is suggested to be dynamically adjusted from 4 Journal of Sensors 0.9 at the beginning to 0.4 at the end of search, and acceleration factors are suggested to be constant as 1.49 or 2. In addition, the particle's population size s can affect the accuracy and efficiency of the algorithm. It should be determined according to the complexity of the fitness function. The execution time of the PSO algorithm increases along with the increase of the population size. In general, the execution time of the PSO algorithm may be longer than that of the L-M algorithm.

Simulation Section
In simulation, the proposed method was applied to track the magnetic target. The parameters of a magnetic target were set as m, θ, φ = 600 Am 2 , 60 deg, 15 deg . The target moved on a path which was parallel to the X-axis, starting at −10 m, 10 m, −1 m and ending at 10 m, 10 m, −1 m , shown in Figure 2. The array consisted of seven sensors. One sensor was located at the origin of the coordinate system, and other sensors were located on the orthogonal axes at 2 m from the origin. The measurement error was set as follows: the average is zero, and the standard deviation is 0.04 nT. In the PSO algorithm, we set c 1 = 2, c 2 = 2, the population size to be 150, and the max iterations to be 200. Then, we used the proposed method to track the magnetic target. The tracking results are shown in Figures 3-7. The positions of the target were calculated by the proposed method and the normal method, shown in Figure 3. It can be seen that the calculated positions were close to the true values. The magnetic anomaly measured by the magnetometers decreased with the cube of the distance, leading to a reduction in the signal-to-noise ratio (SNR). Thus, the location accuracy decreases as the distance increases. Compared to the normal method, the proposed method had better location accuracy when the distance was large. The angles of the magnetic moment were calculated, shown in Figures 4 and 5. It can be seen that the results calculated by the proposed method are better than those by the normal method. The magnitude of the target's magnetic moment is shown in Figure 6. The results show that the calculated values by the proposed method were closer to the true values.

// Initialization
Initialize the position and velocity of each particle with a uniformly distributed random vector, respectively.
where b lo denotes the minimum boundary of searching space or solution space and b up denotes the maximum boundary of searching space or solution space. // Obtain the measured field from the magnetometers for j = 1⋯n do ΔB ij ← read from measured data end for Calculate the fitness of each particle. g x i Initialize the best position of each particle p i ← x i and the best position of the swarm g ← min p i // Main PSO algorithm loop while the criterion is not met do for each particle i = 1⋯s do Update the velocity of each particle: v i ← wv i + c 1 rand 1 p i − x i + c 2 rand 2 g − x i Update the position of each particle: If g x i < p i then Update the particle's best position: p i ← x i end if If g p i < g then Update the particle's best position: g ← p i end if end for end while Output the estimated parameters x t , y t , z t , θ, φ of the mobile target // Evaluate the solution calculated by the optimization method Calculate the quality index QI, and determine whether to accept or reject the solution Output the final parameters x t , y t , z t , θ, φ, m

Journal of Sensors
The quality index was calculated according to (10), shown in Figure 7. It can be seen from Figure 7 that the value of the quality index was small when the solution error of the target was large. Therefore, the quality index can be used to estimate the accuracy of the solution calculated by the dedicated PSO algorithm.
The quality index is used to determine whether the solution is the global optimal solution or the local optimal solution. As we know, the relative error of the global optimal solution is smaller than that of the local optimal solution. In practice, the distance r and the magnitude of the magnetic moment are the major parameters for the tracking system. Thus, the calculated solution can be considered the global optimal solution if the relative error of the distance was less than 15% and the relative error of the moment was less than 50% in the simulation. In this case, we can determine the optimal value of the quality index, which can be used as the threshold value. After many simulations at different points, we determined that the optimal value of the quality index was 0.3.

Experimental Setup.
In order to test the performance of the proposed method in practice, we used it to track a car with a constant velocity which moved along a planned trajectory. The experiment was carried out in Jinshatan Wetland Park in Harbin City, China, where the ambient magnetic    Journal of Sensors activity due to external sources, such as power lines and traffic, was very low. We used the north-east-down (NED) coordinate system in the experiment and located the center of the array at the origin in the experiment. The array consisted of four optically pumped cesium magnetometers (CS-L [Scintrex, Snidercroft Road, Concord, Ontario, Canada]), and each magnetometer was located on the orthogonal axes at 2 m from the origin, as shown in Figure 8. This magnetometer has high sensitivity, and intrinsic noise is about 0.6 pT/√Hz at 1 Hz. In addition, there is also another source of magnetic reading uncertainty in cesium vapour magnetometers called the heading error (±0 2 nT for CS-L). It is a function of the angle of the sensor head with respect to the local magnetic field. In this experiment, we used a static sensor array to track the magnetic target. We did our best to ensure that the magnetometers had the same orientation so that they would have the same heading error. Thus, the heading error was not considered in this experiment. The measured data were acquired at a sampling rate of 2 Hz. An SUV car as the target was driven along the planned trajectory in the horizontal plane. The plan trajectory was parallel to the X-axis with y = 31 05 m. In the experiment, the target started at x = 14 0 m and ended at x = 14 5 m, the Z position of the target was −0.63 m, and the parameters of magnetic moment remained constant at 485 8 Am 2 , 33 8°, 48 7°.

Experimental Results.
Since there were four optically pumped magnetometers in the experiment, we can only obtain three unrelated magnetic anomalies ΔB ij . We can thus estimate only three parameters of the target. Therefore, we set the parameters z t , θ, φ and the magnetic target as known values and applied the proposed method on the real data measured by the magnetometers to estimate the other three parameters x t , y t , m .
The results are shown in Figures 9-11. The estimated positions of the target moving in the X-Y plane are shown in Figure 9, and they are seen to be close to the true values. The localization error at both ends of the trajectory was large compared with that at the middle, because the signal-to-noise ratio (SNR) at the ends was smaller than that at the middle.
After we obtained the optimal solution about parameters x t , y t of the target by the dedicated PSO algorithm, we can use (8) to calculate the magnitude of the magnetic moment. The average magnitude of the magnetic moment of the target was estimated by the magnetic anomaly ΔB ij , as shown in Figure 10. Due to the localization error, the error of the magnetic moment was large. Then, we used the quality index defined in (10) to evaluate the quality of the optimal solution calculated by the proposed method. We can directly determine that the solution was accepted or rejected by the quality index. From (10), the solution is good if the value of the quality index is near 1 and the solution is bad if the value of the quality index is near 0. The threshold value of the quality index for accepting or rejecting the solution should be determined based on the positioning accuracy requirement. The quality index of the solution is shown in Figure 11. We can see that the value of the quality index was small at both ends. It means that the localization error was relatively larger.

Conclusions
In this paper, we presented a method based on a scalar magnetometer array to track a moving magnetic target. In this method, we isolated the magnitude of the magnetic moment from the scalar magnetic anomaly and formulated an optimization problem based on it to determine the position and magnetic parameters of the target. Then, a dedicated particle swarm optimization algorithm was implemented to solve this optimization problem. The PSO algorithm may converge to the local minimum and cannot give the optimal solution. Thus, we define the quality index, which is used as an indicator for determining whether the calculated solution is a global optimal solution. We can improve the accuracy of the tracking of the remote magnetic target based on the quality index. The proposed method was applied to the simulation and the real data to track the magnetic target. The results showed that the tracked trajectory was very close to the true trajectory.

Conflicts of Interest
The authors declare that they have no conflict of interests.  Journal of Sensors