Proceedings of the Institute of Acoustics

 

Robust magnetic dipole trajectory estimation handling false alarms and with noisy measurements

 

Joep de Jong, NATO-STO Centre for Maritime Research and Experimentation (CMRE), La Spezia, Italy
Pietro Stinco, NATO-STO Centre for Maritime Research and Experimentation (CMRE), La Spezia, Italy
Alessandra Tesei, NATO-STO Centre for Maritime Research and Experimentation (CMRE), La Spezia, Italy
Robert Been, NATO-STO Centre for Maritime Research and Experimentation (CMRE), La Spezia, Italy

 

 

1 INTRODUCTION

 

Underwater surveillance with active and passive sonar systems is a well-known problem in the field of underwater acoustics. Moreover, recent advances in sensor technology open the possibility of integration of new generation, low Size, Weight, Power and Cost (SWaP-C), quantum magnetic sensors into underwater nodes for area surveillance through magnetic anomaly detection. In recent years, CMRE developed a set of bottom stations equipped with passive acoustic sensors, such as acoustic vector sensors and compact volumetric arrays, as well as total-field quantum magnetic sensors. Data and information fusion of acoustic and magnetic data are exploited to increase the detection capability of a single bottom node as well the localization accuracy at the Command and Control (C2) center where acoustics and magnetic detections are collected.

 

This work establishes a method for estimating the trajectory of a moving magnetic dipole using a network of bottom nodes equipped with total-field magnetometers. With the limited communication capabilities in the underwater domain, the sharing of raw total-field signals between the bottom nodes for noise removal using multi-channel coherence 7 is prohibited, consequently, each node has to decide, locally, whether a dipole is detected and communicate the bare minimum of information to the central processing unit to estimate the trajectory of the target. The research presented in this paper will focus on a robust algorithm able to estimate the target trajectory of the moving dipole with noisy measurements of the dipoles detected and in the presence of false alarms. The proposed algorithm is based on a fitness function, which allows more freedom in finding a trajectory and assigns a score to a trajectory that indicates, for each sensor, how the estimated trajectory aligns with the measured values of the dipoles at each node.

 

2 BACKGROUND

 

At sufficiently large distances, the magnetic field of a moving ferromagnetic object can be modeled as a magnetic dipole. The magnetic field of a dipole is given by

 

 

where B is the magnetic field in T, m the dipole moment in Am 2, r the position vector of the dipole and the sensor in m, and µ 0 the permeability of free space in N / A 2. An anomaly caused by the relative movement between a dipole and a magnetometer can be characterized by the Anderson functions 9: a basis for the magnetic field of a dipole moving in a straight line with constant velocity. They are defined as

 

 

where θ is a dimensionless parameter indicating the position on the line. The representation of the magnetic field of a dipole on a line in terms of three Anderson functions is given by

 

 

where C_n are the coefficients of the Anderson functions:

 

 

A single magnetometer can detect a moving dipole target, by looking at deviations from a reference field^1,4. As the magnetic field of a dipole decays with a third power of the distance, a dip (or peak) in the magnetic field occurs when the target is at its closest point of approach (CPA) at time t. This anomaly can be characterized by the Anderson functions (2). The shape of the anomaly determines the characteristic time τ : the ratio between the (unknown) distance | r | to the CPA and the (unknown) velocity v of the target:

 

 

The Anderson functions are used as a matched filter to detect the anomaly in the measured field. A multi-scale approach 5 is used to detect anomalies of different shapes (different distances and velocities), leading to a search space of possible time shifts T = [ t ₁ , t ₂ , . . . , t _n ] in s and characteristic times T = [ t ₁ , t ₂ , . . . , t _m ] in s . Related work describes this process in detail^9,8. The pseudo-energy of the magnetic signal E can be defined as the area under the curve of normalized Anderson functions, given by the sum of the Anderson coefficients. Maximizing the pseudo-energy over the space T × T gives an estimate of the CPA (6). The closest point of approach (CPA) obtained from a single sensor i can therefore be expressed as a triplet (t_i, T_i, E_i).

 

 

 

3 PROBLEM STATEMENT

 

The problem to be addressed is the estimation of the trajectory of a dipole using data from multiple scalar magnetometers. The assumption is made that a dipole moves in a straight line (heading ϕ ) with constant veloc- ity v . We consider bottom nodes in a surveillance area with a given bathymetry and hence at different depths. Given n magnetometers, positioned at S_i ∈ Ω ⊂ R³ for i = 1 . . . n , where Ω is the domain of interest. The measurements from the n magnetometers are:

 

t = ( t ₁ , . . . , t _n ) , τ = ( t ₁ , . . . , t _n ) and E = ( E ₁ , . . . , E _n ) .

 

The problem can be reformulated as a geometric optimization problem with input time t _i and characteristic time τ _i . Note that the relation between the characteristic time and the distance to the closest point of approach is given by r _i = τ _i v . Therefore, under noise-free conditions, the CPA lays on a sphere with radius τ _i v around sensor i , and the direction of movement ϕ is tangent to this sphere. A unique solution for this tangent line exists, given measurements of at least three sensors 2. The goal is to find, in the presence of noise and false alarms, a pair of ( v, ϕ ) for which (i) a tangent line to all spheres with radius τ_iv exists, and (ii) the times t i are consistent with the time it takes to reach the different CPAs, as illustrated in figure 1.

 

 

Figure 1: Trajectory of a target being tangent to circles with radius r_i = τ_i v.

 

4 ALGORITHM

 

In the presence of noise, the existence of a tangent line to all spheres with radius τ_iv is not guaranteed. To address this, we aim for the line most tangent to the spheres, achieved through a fitness function that scores trajectories based on the parameters v, ϕ, γ₀, h, ϵ_τ and ϵ_t, where v in m s⁻¹ and ϕ in ^◦ are the velocity and the heading of the target, γ₀ is the sign of the north/south ambiguity for the global closest point of approach x₀ (GCPA), h in m is the altitude of the dipole, and ϵ_τ and ϵ_t in s are corrections in respectively the characteristic time and the time, used to find the GCPA. These parameters encompass all degrees of freedom in the problem: position x₀, movement direction ϕ, velocity v of the target and time. 

 

The algorithm starts with a GCPA (t₀, τ₀), determined initially. One straightforward method is to select the sensor S₀ with the smallest τ value, located at (x₀, y₀, z₀). Then, a radius r₀ is chosen, and movement occurs orthogonally to ϕ, considering the north/south ambiguity γ₀. The altitude of GCPA is set to h, using ψ_i = arcsin((h − z_i)/r_i) to find the projection of r_i  on the 2D plane:

 

 

Assuming measurements (t_i, τ_i, E_i) are sorted on their characteristic time, τ₁ ≤ τ₂ ≤ . . . ≤ τ_n, the GCPA is defined as the corrected version of the measurements with the smallest τ : (t₀, τ₀) = (t₁ + ϵ_t, τ₁ + ϵ_τ ), where ϵ_t and ϵ_τ correct the matched filter’s estimation of τ₁ and t₁, to ensure that moving with γ₀r₀ gives a point on the line. The position of the GCPA C^γ0

 

 

and the trajectory of the dipole is given by

 

 

The sensor positions S_i projected on the trajectory (9) yield the distances d_i ^γ0 along the movement line. The projection in orthogonal direction gives the projected distance from the sensor to the CPA r̂ ^′_i^γ0:

 

 

 

Dividing this distance by the velocity gives the estimated time t̂_i ^γ0 and the position of the CPA C_i ^γ0 estimate of r_i  can be obtained by combining r̂ ^′_i^γ0 with the vertical displacement:

 

 

4.1 Fitness Function

 

Two fitness functions can be defined. The first represents the radial error expressed in meters m:

 

 

where r_i = τ_iv is measured distance and r̂ _i^γ0  is the estimated distance. Different from related work⁶, we do not rely on the maximum total-field values, and a L1-norm, instead of L2, is used to prevent outliers from dominating the radial error. The second fitness function represents the tangent error in seconds s:

 

 

The fitness functions assign a score to a trajectory, given the parameters v, ϕ, γ₀, h, ϵ_τ and ϵ_t. The values range between 0 and ∞, where 0, the best score, is attained when the measured t_i and τ_i are consistent with the estimated t̂_i and T̂_i for all sensors and the given parameters. In a noise-free scenario without false alarms, there is a unique solution for the parameters v and ϕ that minimizes the fitness functions. In the presence of noise and false alarms, a zero fitting error is not guaranteed. The goal is to find the best estimate for (v, ϕ, γ₀), h, ϵ_τ and ϵ_t that minimizes the fitness functions. 

 

The problem is described using two fitness functions and is therefore multi-objective. However, it can be made single-objective by expressing the fitness functions in the same unit. This can be done by either dividing the radial error J_ri (m) or multiplying the tangent error J_ti (s) by the velocity v (m s⁻¹). Since the errors represented are orthogonal, the total error J_i can be expressed as the Euclidean distance between the two errors. The global fitness function J is defined as the sum of the fitness functions for all sensors.

 

 

4.2 False Alarm Classification

 

An algorithm has been developed to improve accuracy by implementing a simple false alarm classifier. This classifier operates without prior knowledge of the target’s movement direction, relying on the principle that the Closest Point of Approach (CPA) time should align with the target’s velocity. 

 

Initially, a minimum velocity threshold is set. The classifier evaluates whether a target takes an excessively long time to reach the CPA between two sensors. This time depends on the target’s heading ϕ, with an upper bound determined by the distance between sensors divided by the minimum velocity. Input parameters moving too slowly are removed from the data. Then, the signal energy is filtered using a Robust z-score method³, which involves computing a ratio tied to the magnetic moment. This ratio, independent of velocity, is derived by multiplying signal energy (E_i) by the cube of the characteristic time (τ³). A threshold of 3 is set to filter out false alarms. 

 

While effective, the classifier’s performance can be significantly enhanced with prior estimates of the heading ϕ: stricter values of the minimum and maximum times between CPAs improve the filter. 

 

4.3 Optimizing Fitness

 

The fitness function J_i as defined in (13) represents the error that is made in the estimated CPAs, and the measured CPAs. To minimize the fitness function for the unknowns ϕ, v, h, γ₀, ϵ_τ and ϵ_t, a gradient descent algorithm is applied. For this, the partial derivatives are calculated in closed form and are given in Appendix A. The algorithm runs in a 5-dimensional space, where the unknowns are the parameters v, ϕ, h, ϵ_τ and ϵ_t. It is suggested to run the algorithm in parallel for both γ₀ = 1 and γ₀ = −1. Once either of the algorithms converges, the best estimate is returned, otherwise, we limit the number of iterations to prevent the algorithm from running indefinitely. It is guaranteed that the algorithm converges to a local minimum. To have a higher chance of finding the global minimum, the algorithm can be run multiple times with different starting points. Initialization can be done with a guess if available, otherwise, the fitness function can be evaluated for a range of ϕ and v to find a good starting point.

 

5 SIMULATION

 

For the simulation, a domain of 5000m × 5000m with n = 3 . . . 9 sensors is used. The sensors are positioned at random locations within the domain, or on a circle with a radius of 1500m forming a polygonal shape. The sensors are positioned at a depth of 100m. Figure 2 shows some of the different sensor distributions with 7 sensors, used in the simulations. In practice, the input parameters, consisting of the

 

 

Figure 2: Examples of different sensor distributions used in the simulations.

 

CPAs of each sensor (t_i, τ_i, E_i) are obtained from the measured signal after applying the matched filter, as described in Section 2. Under low SNR conditions, the main uncertainty in the estimation of the CPAs is in the characteristic time τ . The time t is less prone to errors, as the matched filter is more accurate in estimating the time of the CPA. In the simulation, we take the true values of t and τ and add noise to simulate measurement error. The noise is sampled from a uniform distribution at different noise levels. 

 

To generate false alarms, a sensor k is randomly selected and a time t_k is randomly selected within a given time window. The characteristic time τ_k is replaced with a random value within the range of existing τ values and the energy is randomly set to a value within the range of the energy values. 

 

 

Figure 3: Real and estimated trajectory us- ing a 10-sensor setup with 2 false alarms.   

                                                                                            

 

Figure 4: The effect of false alarmson the estimation of the heading ϕ .

 

6 RESULTS

 

The outcome ideally involves the trajectory being tangent to all spheres of radius τ_iv around sensor i, excluding those identified as false alarms. This objective is equivalent to minimizing J_i for all non-false alarm sensors. Figure 3 illustrates the results of a simulation conducted with ten magnetometers, two of which yielded false alarms. The algorithm has successfully identified a solution that aligns well with accurate measurements but performs poorly on the false alarms (9 and 10). In this solution, the estimated trajectory tangentially intersects all circles with radii r_i = vτ_i for sensors without false alarms, closely approximating the true trajectory. 

 

To assess the accuracy of parameter estimation for ϕ, we apply thresholding to the absolute deviation. Specifically, a threshold of 5^◦is set for the heading ϕ. Figure 4 displays the outcomes of simulations with varying numbers of sensors and false alarms. Each simulation is ran 1000 times. The results indicate the algorithm’s robustness against false alarms, with a high percentage of experiments falling within the threshold. Moreover, the algorithm’s performance tends to improve with an increase in the number of sensors. Another metric used to evaluate the algorithm is the RMSE of the estimated CPAs. The CPAs are calculated using the estimated parameters v, ϕ, γ, h, ϵ_τ and ϵ_t and the measured times t_i and distances τ_i

 

 

Figure 5: Proportion of samples with an RMSE below a certain value.

 

Figure 5a shows what proportion of 1000 runs has an RMSE of the estimated CPAs less than a certain value. A setup with seven sensors in a polygonal setup and different numbers of false alarms is used. As the count of false alarms increases, the RMSE tends to rise. A priori filtering of the false alarms with the method presented can significantly improve the performance of the algorithm. With either zero or one false alarm, the majority of experiments exhibit an RMSE below 50m. Moreover, around 85% of the experiments achieve an RMSE below 200m with seven sensors and two false alarms.

 

The performance of the algorithm under various noise levels is depicted in Figure 5b, employing eight sensors positioned randomly and with no false alarms. As the noise levels increase, there is an increase in the RMSE. It is important to highlight that, due to the random distribution of sensors, certain simulations approximate worst-case scenarios, with poor sensitivity to the fitting metrics. This aspect explains the high RMSE observed in certain simulations. 

 

Under false alarms and increasing noise levels, the following observations are made: In the absence of noise, with two false alarms present, approximately 93% of the simulations adhere to the 5-degree threshold (Figure 4). With a 10% noise level in τ and 10-second noise in t, this percentage decreases to 90%. 

 

7 CONCLUSION

 

We presented a method to estimate the trajectory of a dipole moving in a straight line using a network of nodes equipped with total-field magnetometers. The described procedure requires minimum communication between the underwater nodes: no raw data is required to be exchanged: instead only a few bytes containing the locally computed information about CPA, consisting of the time t_i , the characteristic time τ_i  and the energy E_i , must be communicated to the C2 data fusion nodes. 

 

The central C2 data fusion node can estimate the trajectory of the target by solving a geometric problem. The trajectory must be tangent to the line between the sensor and its CPA. The length of this ”radial” line is approximately given by r_i = vτ_i . Two fitness functions are presented: the radial error, defining the error in the length of the radial line, and the ”tangent” error, defining the error between the measured time t_i and the estimated time t̂_i^γ0 at which the CPA occurs. A global fitness function is presented and a gradient descent algorithm is applied to it by using the partial derivatives in closed form. Minimization of the fitness value determines the trajectory with minimal radial- and tangential error. 

 

Simulations show that the algorithm performs very well under false alarms. The algorithm is also robust to significant noise levels in the measurements. It is also shown that increasing the number of sensors improves the performance of the algorithm. 

 

Future research will focus on data fusion between acoustic and magnetic data to increase the probability of detection of a bottom node as well the localization accuracy achievable at a Command and Control (C2) station where acoustic and magnetic detections are collected to estimate the trajectories of the detected targets.

 

A APPENDIX

 

Partial derivatives of the Radial Error

 

 

Partial derivatives of the Tangent Error

 

 

REFERENCES

 

1. P. Alken et al. International Geomagnetic Reference Field: The thirteenth generation. Earth, Planets and Space, 73(1):1–25, December 2021.

2. C. P. Du, M. Y. Xia, S. X. Huang, Z. H. Xu, X. Peng, and H. Guo. Detection of a moving magnetic dipole target using multiple scalar magnetometers. IEEE Geoscience and Remote Sensing Letters, 14(7):1166–1170, 2017.

3. E. Elamir. Mean absolute deviation about median as a tool of explanatory data analysis. In Proceedings of the World Congress on Engineering, volume 1, 2012.

4. B. Maus et al. Emag2: A 2–arc min resolution Earth magnetic anomaly grid compiled from satellite, airborne, and marine magnetic measurements. Geochemistry, Geophysics, Geosystems, 10(8), 2009.

5. Boris Ginzburg, Lev Frumkis, and Ben-Zion Kaplan. Processing of magnetic scalar gradiometer signals using orthonormalized functions. Sensors and Actuators A: Physical, 102(1):67–75, 2002.

6. Richard J. Kozick and Brian M. Sadler. Algorithms for tracking with an array of magnetic sensors. In 2008 5th IEEE Sensor Array and Multichannel Signal Processing Workshop, pages 423–427, 2008.

7. C. E. Lucas and R. Otnes. Noise removal using multi-channel coherence. Defence R&D Canada–Atlantic, 2010.

8. Arie Sheinker, Boris Ginzburg, Nizan Salomonski, Phineas A. Dickstein, Lev Frumkis, and Ben-Zion Kaplan. Magnetic anomaly detection using high-order crossing method. IEEE Transactions on Geoscience and Remote Sensing, 50(4):1095–1103, 2011.

9. W. Michael Wynn. Detection, localization, and characterization of static magnetic-dipole sources. In Detection and Identification of Visually Obscured Targets, pages 337–374. Routledge, 2019.