A Robust Tracking Method for Multiple Moving Targets Based on Equivalent Magnetic Force

A ferromagnetic vehicle, such as a submarine, magnetized by the Earth’s magnetic field produces a magnetic anomaly field, and the tracking of moving targets can be realized through real-time analysis of magnetic data. At present, there are few tracking methods based on magnetic field vectors and their gradient tensor. In this paper, the magnetic field vector and its gradient tensor are used to calculate equivalent magnetic force. It shows the direction of the vector between the detector and the tracking targets for controlling the direction of motion of the detector and achieving the purpose of tracking. Compared with existing positioning methods, the proposed method is relatively less affected by instrument resolution and noise and maintains robustness when the velocity vectors of multiple magnetic targets change randomly.


Introduction
A magnetometer in geophysical applications is a precision instrument used to observe magnetic anomalies caused by rocks, ores, and ferromagnetic objects (such as submarines) which are magnetized by the geomagnetic field [1][2][3]. It is the basis for studying the geological structure, mineral resource distribution, and magnetic characteristics of ferromagnetic objects. With the development of modern physics, many institutions have developed magnetometers based on proton precession, fluxgate, optical pump, and superconductivity [4][5][6][7], and the accuracy of magnetic measurement has continuously improved. According to physical quantities to measure, magnetometers can be divided into three categories: total field magnetometers for the magnitude of the magnetic field, vector magnetometers for three components of the magnetic field [8,9], and tensor magnetometers for magnetic gradient tensor [10,11]. Among them, the magnetic gradient tensor provides rich information and suppresses time-domain interference of the geomagnetic field. It is especially suitable for measurement on mobile platforms and has become an international research frontier [3,10].
An essential application of magnetometers is the real-time tracking of magnetic targets [12]. The current development trend mainly has two aspects: firstly, the magnetic targets to be tracked are more and more abundant, and their motion states are more and more complex [13,14]; secondly, the tracking methods have gradually developed from the initial use of the total field magnetometer to the combination of multiple magnetometers [15,16], such as the combination of vector magnetometer and tensor magnetometer. This is because the total field magnetometer can only recognize the targets' existence within its detection range through changing magnetic field value from the magnetic anomaly, but it cannot track moving targets due to lack of directivity [17,18]. For example, the responsive search submarine work is to establish a specific search model based on the total field magnetometer or other instruments in a special search scene, and after determining the search range, improve the submersible search efficiency by studying different routes of the detector Micromachines 2022, 13, 2018 2 of 12 or other environment-affected parameters [12,[19][20][21][22][23]. The combined tracking method of vector magnetometer and tensor magnetometer mainly adopts Euler deconvolution [15,24]. The traditional Euler deconvolution method mainly establishes a linear equation system about the position vector through the potential field relationship between the magnetic field and its gradient tensor. However, there is a problem with numerical calculation stability in tensor inversion. Nara et al. (2014) have pointed out that when the magnetic moment direction is perpendicular to the position vector, the magnetic gradient tensor matrix becomes a singular matrix [25]. This situation is more likely to occur where multiple magnetic targets coexist, with more errors involved in solving results [26]. To better ensure the stability of inversion, Nara et al. have proposed the Euler deconvolution method based on generalized inversion, which has partially solved the numerical stability problem [25]. Yin et al. (2020) have re-analyzed the relationship among position vector, magnetic field vector, and magnetic gradient tensor to avoid inversion calculation of the magnetic gradient tensor [27], but it is still affected by magnetic field noise.
The goal of this paper is the real-time tracking of moving targets in multi-target fields. To that end, we offer an analysis method based on the equivalent magnetic force that can successfully track moving targets in such fields. The tracking performance is inevitably impacted by various noises during motion tracking, particularly the inaccuracy of magnetic field measurement while the detector is moving. The equivalent magnetic force method is more advantageous than the currently employed positioning methods. Because of its high robustness, ease of use, and calculation stability, it can effectively weaken the impact of various noise and measurement errors in the tracking process and guarantee the real-time performance of motion tracking. Additionally, where multiple magnetic targets exist simultaneously, this method can accurately track the moving target with the highest equivalent magnetic force.

Introduction to Tracking Models
When a magnetic target and a detector are both moving continuously, motion tracking is the process of detecting and analyzing real-time changes in the target's magnetic field and immediately adjusting the detector's motion direction. Typically, the magnetic detector keeps a large distance from the tracked target and may not be in the same motion plane. This paper uses a submarine with magnetic properties as an example for detection and tracking to simplify the motion tracking model, as shown in Figure 1. In the tracking process, it can be considered that the target's initial position is randomly distributed within a specific range. The target has a maximum speed max and its initial speed follows a uniform distribution on [0 , max ]. The heading of the target fol- We equivalent the motion of the magnetic target into the following four models: Model 1: The heading and the speed remain constant during the movement of the mag- It is possible to generate independent random motion trajectories by positioning multiple moving targets in the same plane and allowing them to move independently. The detector performs motion tracking in a plane higher than the moving target, so the target is equivalent to a moving magnetic dipole [28]. The measured magnetic field vector and magnetic gradient tensor data change in real-time due to the changing magnetic field as the target moves. The orientation information of the magnetic target and the direction vector of the detector at the next moment can be obtained by analyzing and processing the detection data at the current moment. The detector adjusts its motion direction to the direction vector toward the target.
Practically speaking, the magnetic anomaly is small in comparison to the total geomagnetic field; therefore, the motion platform's magnetic field vector is easily influenced by geomagnetic field fluctuations, leading to data inaccuracy. Additionally, the noise of the moving platform and other environmental interference makes tracking more difficult and affects the accuracy of computation results. As a result, a data analysis method that is less influenced by error interference and has high robustness is required.
In the tracking process, it can be considered that the target's initial position is randomly distributed within a specific range. The target has a maximum speed v max and its initial speed follows a uniform distribution on [0, v max ]. The heading θ of the target follows a uniform distribution on − η 2 , η 2 (where η 2π) [17]. We equivalent the motion of the magnetic target into the following four models: Model 1: The heading and the speed remain constant during the movement of the magnetic target. Model 2: The heading remains constant, and the speed changes randomly during the movement of the magnetic target. Model 3: The heading changes randomly, and the speed remains constant during the movement of the magnetic target. Model 4: The heading and the speed change randomly during the movement of the magnetic target. Convert the coordinate system to a Cartesian coordinate system and keep it constant throughout the detection process. Assume the speed of the detector to be v max , which is the maximum speed of the target, and the sampling interval to be t 0 . The direction vector s of the detector at the next moment is (L, M, N) in three dimensions. The detector's displacements in x and y directions in an interval t 0 correspond to ∆X and ∆Y, respectively, shown in Equations (1) and (2): The detector tracks the moving target in the XOY plane. Assume the initial position of the detector to be (X 0 , Y 0 , Z 0 ) and the magnetic target to be (x 0 , y 0 , z 0 ). Then at the next moment t 1 , the former will become (X 1 , Y 1 , Z 1 ) = (X 0 + ∆X, Y 0 + ∆Y, Z), and the latter will become (x 1 , y 1 , z 1 ). At the same time, the relative position relationship between the detector and the moving target has changed. At the next moment t 2 , the direction vector s of the detector will be recalculated, and the above process will be repeated. Finally, the detector realizes the motion-tracking function in the iterative calculation process.

Introduction to Tracking Methods
From Ampere's Law, the elementary force on a current element IdI in the presence of a magnetic field B is given by dF = I(dI × B) [29]. We can further organize the formula to obtain Equation (3) [30]: where M r is the remanent magnetic moment, µ 0 is the magnetic permeability, χ is the magnetic susceptibility of the test mass, and V is the volume of the test mass.
When in the presence of a magnetic force sensor, Equation (3) can be used to calculate the precise magnetic force of the magnetic material. Otherwise, when only the magnetic field vector and the magnetic gradient tensor are known, and the target's magnetic properties are unknown, Equation (4) can be used to calculate the equivalent magnetic force instead.
where B = B x ; B y ; B z is the magnetic field vector, the vector F = F x ; F y ; F z represents the equivalent magnetic force, and Equations (5) and (6) show its components in the x direction and the y direction, respectively.
The magnetic field formed by a magnetic dipole at a single observation point can be expressed by Equation (7).
where M is the magnetic moment of the magnetic dipole, r is the source-to-sensor position vector, µ is the air magnetic permeability, i, j, k = 1, 2, 3 represent x, y, and z in the Cartesian coordinate system. The magnetic gradient tensor G is defined as the gradient of the magnetic field vector B [31], which has nine components and can be represented by the following matrix.
The magnetic gradient tensor G can be expressed by the scaled moments and direction tensor [32], each component of which can be represented by Equation (9). Similarly, each component of the magnetic field vector B in Equation (7) can be represented by Equation (10).
where M B i = m i r 3 and M G i = 3m i r 4 are the components of M B and M G respectively, M B represents the scaled moments of B, and M G represents the scaled moments of G. N B ij = 3n i n j -δ ij and N G ijk = 5n i n j n k − δ ki n j + δ kj n i + δ ij n k are the components of the direction tensors N B and N G , respectively, where N B is a second-order tensor and N G is a third-order tensor. Both N B and N G are functions of the direction cosine (n 1 , n 2 , n 3 ) of the position vector r.

of 12
Substitute Equations (9) and (10) into Equation (4) and derive the equivalent magnetic force, which is the product of the scaled moments and the orientation tensor as Equation (11).
In Equation (11), the principle of tensor analysis [33] is applied to express the equivalent magnetic force as the product of tensors. Then the direction tensor is contracted once to obtain the new direction tensor S of the equivalent magnetic force F. As shown in Equation (12), S is also a third-order tensor.
klm g k g l g m = N G ij N B klm g i g k g m = S ikm g i g k g m (12) In this paper, the tracking process is discussed in a Cartesian coordinate system with both base vectors normalized and orthogonal. So when contracting the tensor, it is not necessary to distinguish between the upper and lower indexes of base vectors and the obtained dummy index, so each component S ijk can be expressed as Equation (13): The direction tensor S represents the direction of the position vector r. However, due to the objective existence of the magnetic moment direction, the calculated direction vector s F from the equivalent magnetic force method is influenced by M B , M G , and S at the same time. The direction vector s F of the detector is obtained after a little deviation, but there is a real-time relative motion relationship between the detector and the moving magnetic source during the motion tracking process, and there is no cumulative error, so it has little impact on the overall tracking effect.
During the tracking process, the detector will form a tracking trajectory following the moving target. By drawing error bands of the tracking trajectory, the effect of the magnetic moment direction on the direction vector s F is clearly demonstrated. The process and principle of drawing error bands are shown in Figure 2.

Motion Tracking Simulation
Each moving target is considered a separate magnetic submarine, which generates The relative position relationship at a certain moment is used to explain the drawing principle. The blue trace depicts all the positions that the detector might reach at the next moment after being affected by all possible magnetic moment directions. Among these tracking points, respectively find the points that are farthest from the actual pointing (red arrow) vertical distance from the detector to the tracking target. Both sides are the upper and lower error band sample points. As shown in Figure 2, the black dots are the positions of the detector, and the red dots are the positions of the tracking target. It can be observed that the tracking trajectories generated by the equivalent magnetic force analysis method are always within the drawn error band.

Motion Tracking Simulation
Each moving target is considered a separate magnetic submarine, which generates random trajectories according to the motion models.

Motion Tracking Simulation
Each moving target is considered a separate magnetic submarine, which generates random trajectories according to the motion models.   Among all the moving targets, the direction vector points to the one with the greatest equivalent magnetic force or the closest relative distance from the detector, which is thus tracked in the subsequent process. We set four moving targets with the same magnetic field intensity, as shown in Figure 3. In the initial position, target 3 has both the closest relative distance and the greatest equivalent magnetic force compared with others. Therefore, at the early stage of tracking, the detector will move towards target 3 according to the direction vector s F . However, as the tracking progresses, the equivalent magnetic force of target 4 gradually becomes the greatest, and at a specific time, the next target changes into target 4. The motion trajectories of different targets can occasionally cross because they move in completely distinct ways. Since the detector always tracks the moving target with the highest equivalent magnetic force, the tracking target is not lost. As shown in Figure 3, when the motion trajectories of target 2 and target 4 crossed, the tracking target of the detector switched from target 4 to target 2. Nara et al. (2006) have demonstrated the principal formula of Euler's deconvolution [15].

Comparative Analysis of Noise Interference
In the calculation process, the inverse of the magnetic gradient tensor G needs to be solved.
When the magnetic moment vector M is perpendicular to the position vector r, the magnetic gradient tensor matrix changes into a singular matrix, and the condition number approaches infinity, according to a discovery by Nara et al. [25]. Considering this situation, they proposed the Euler deconvolution method based on Moore-Penrose generalized inverse. When the magnetic gradient tensor matrix is singular, the generalized inverse method can be used to obtain the unique solution of the position vector.
According to Zhang's theory [34], the condition number is a key index to quantify numerical stability. During motion tracking, the angle between the magnetic moment vector M and the position vector r inevitably approaches 90 degrees in some cases. Since the magnetic gradient tensor matrix G is not singular then, corrections cannot be achieved using the generalized inverse method. In such cases, however, the condition number is relatively large, and the equation for the position vector is pathological, which means that when the initial value of the magnetic field vector B is slightly disturbed, the solution of the equation changes significantly.
Based on these works, Yin et al. focused on the relationships between the analytical formula and the gradient tensor of the magnetic field vector to avoid the inversion problem associated with traditional Euler deconvolution methods [27]. The non-inverting Euler method proposed by Yin et al. can be organized as Equation (15).
where λ = (λ min , λ med , λ max ) represents the eigenvalues of the magnetic gradient tensor matrix G and E is a 3 × 3 unit matrix. Compared with Equation (4), which represents the equivalent magnetic force method, the Equation (15) still introduces magnetic disturbances through its second term, which further affects the tracking performance in actual situations.
The actual tracking process generates a variety of complex motion trajectories. For different methods, their tracking performances, especially in a superimposed noise environment, show a contrast in robustness. The motion tracking time t is set to 100 s, and a complex motion trajectory is obtained in the simulation. The magnetic field vector B and its gradient tensor G are then superimposed with a Gaussian white noise having a signal-tonoise ratio of 4 and 10, respectively. As seen in Figure 4, when the same noise interferences are superimposed, the tracking performances of the Euler deconvolution method and the generalized inverse method are both significantly impacted, while the equivalent magnetic force method and the non-inverting Euler method perform much better.
The numerical calculation process becomes unstable and the tracking effect is not satisfactory since the first two methods require inversion of the magnetic gradient tensor matrix G. In contrast, the latter two methods do not have this problem. However, a careful comparison of these two methods will discover that the robustness of the non-inverting Euler method is still significantly lower than that of the equivalent magnetic force method, as shown in Figure 5.
Keep the noise level and run the simulation 100 times at random. Then the noninverting Euler method and the equivalent magnetic force method may each provide 100 tracking trajectories. Figure 5 plots all tracking trajectories collectively. By comparison, it is evident that the equivalent magnetic force method is more stable and reliable than the non-inverting Euler method, which further verifies the high robustness of the proposed method in this paper.
to-noise ratio of 4 and 10, respectively. As seen in Figure 4, when the same noise interferences are superimposed, the tracking performances of the Euler deconvolution method and the generalized inverse method are both significantly impacted, while the equivalent magnetic force method and the non-inverting Euler method perform much better. The numerical calculation process becomes unstable and the tracking effect is not satisfactory since the first two methods require inversion of the magnetic gradient tensor matrix G. In contrast, the latter two methods do not have this problem. However, a careful comparison of these two methods will discover that the robustness of the non-inverting Euler method is still significantly lower than that of the equivalent magnetic force method, as shown in Figure 5.
Keep the noise level and run the simulation 100 times at random. Then the non-inverting Euler method and the equivalent magnetic force method may each provide 100 tracking trajectories. Figure 5 plots all tracking trajectories collectively. By comparison, it is evident that the equivalent magnetic force method is more stable and reliable than the non-inverting Euler method, which further verifies the high robustness of the proposed method in this paper.
The mean square error (MSE) between the motion trajectory and the tracking trajectory is used as an index to evaluate the tracking performance. In order to prove that the equivalent magnetic method has a better tracking performance under any complex trajectory, the non-inverting Euler method is still used as the comparison method. Following Model 4, 100 arbitrary motion trajectories are randomly generated, and the tracking process is repeated 100 times at the same noise level, calculating the MSE each time. As shown in Figure 6, the mean value of the MSE for any complex motion trajectory using the noninverting Euler method is larger than the equivalent magnetic force method.  The mean square error (MSE) between the motion trajectory and the tracking trajectory is used as an index to evaluate the tracking performance. In order to prove that the equivalent magnetic method has a better tracking performance under any complex trajectory, the non-inverting Euler method is still used as the comparison method. Following Model 4, 100 arbitrary motion trajectories are randomly generated, and the tracking process is repeated 100 times at the same noise level, calculating the MSE each time. As shown in Figure 6, the mean value of the MSE for any complex motion trajectory using the non-inverting Euler method is larger than the equivalent magnetic force method.

Analysis of Multi-Source Coupling Interference
When multiple magnetic targets coexist in the detection area, the overall magnetic field is the superposition of the magnetic fields from all the targets. Assuming that there are two magnetic targets in the detection area, the measured magnetic gradient tensor G and magnetic field vector B are, respectively, expressed as Equations (16) and (17) according to the superposition principle. Suppose the eigenvalue of target 1 is 1 = ( min 1 , 1 , max 1 ) , the eigenvalue of target 2 is 2 = ( min 2 , 2 , max 2 ) , and the eigenvalue of Equation (16)

Analysis of Multi-Source Coupling Interference
When multiple magnetic targets coexist in the detection area, the overall magnetic field is the superposition of the magnetic fields from all the targets. Assuming that there are two magnetic targets in the detection area, the measured magnetic gradient tensor G and magnetic field vector B are, respectively, expressed as Equations (16) and (17) according to the superposition principle. Suppose the eigenvalue of target 1 is λ 1 = λ 1 min , λ 1 med , λ 1 max , the eigenvalue of target 2 is λ 2 = λ 2 min , λ 2 med , λ 2 max , and the eigenvalue of Equation (16) is λ s = λ s min , λ s med , λ s max . It is unreasonable to consider λ s = λ 1 + λ 2 .
The numerical calculation formula of the Euler deconvolution method can be written as Equation (18): Each term in Equation (18) will change accordingly due to the coexistence of multiple magnetic targets, assuming that Substitute Equation (19) into Equation (18), and then we can obtain Equation (20).
In the field where two magnetic targets coexist, according to the superposition principle, each component needs to be adjusted, for example, G 11 is as shown in Equation (21) Here we analyze a small polynomial that appears during Euler deconvolution, such as G 11 B x : Each item of ∆ r is an interference item caused by the interaction between two targets. The impact will be stronger when there coexist four magnetic targets.
The calculation formula for the equivalent magnetic force method can be written as Equation (24).
Similar to Equation (20), Equation (24) can also be written as follows: We also analyze a small polynomial in the equivalent force method, such as G xx B x : The value in ∆ f corresponds to the interference term caused by the interaction between two magnetic targets. A comparison of Equations (23) and (27) shows that although the equivalent magnetic force method also has interference terms, they are much less than those of the Euler deconvolution method. Note that the influence of interference in |G| has not been considered in this comparison.
Additionally, in multi-target situations, the second term of the non-inverting Euler method is more strongly impacted by the superposition of various error disturbances in the targets' magnetic field vectors, as shown in Equation (28). Additionally, λ s = λ 1 + λ 2 is considered unreasonable when using Equation (15) for calculation. Therefore, the noninverting Euler method is incorrect when multiple magnetic targets coexist.
Theoretical analysis can further demonstrate that the equivalent magnetic force method retains excellent robustness in multi-target situations.

Discussion and Conclusions
The proposed method in this paper can be used to track magnetic targets more accurately after the total field magnetometer determines the existence of magnetic anomalies in the detection range. Remember that if the vertical distance between the detector and the target is set too far apart, the direction of the equivalent magnetic force s F may be primarily biased toward the z direction and the direction vector s F may be more easily affected and become inaccurate. This is because both the magnetic moment direction M and the direction tensor S have impacts on the direction vector s F . For instance, when the magnetic inclination of the detection area is relatively small, the vertical distance above should be reduced to ensure the accuracy of the direction vector s F in the XOY plane because the magnetic moment direction has a greater impact on the x direction and y direction relative to the z direction.
In addition, force balancing issues may occur when the detector is affected by multiple targets. However, on the one hand, there are many disturbances in the natural detection environment, so a constant balance is unlikely to exist. On the other hand, due to the influence of the magnetic moment direction vector M, the direction vector s F will deviate from the equilibrium position, further ensuring the robustness of the proposed method.
The method in this paper is based on the direction vector of the equivalent magnetic force, and the numerical solution process is stable. It significantly guarantees high robustness and reduces the effects of noise disturbances, particularly the measurement inaccuracy of the magnetic field vector, on the tracking performance. In future research, we will analyze the numerical relationship of equivalent magnetic forces more accurately and deeply to obtain more information about moving magnetic targets, such as their magnetic magnitude.

Patents
This section is not mandatory but may be added if there are patents resulting from the work reported in this manuscript.  Data Availability Statement: This study does not report any data.